7 


Marine  Physical  Laboratory 


AD-A270  570 


UMPE: 


The  University  of  Miami  Parabolic  Equation  Model 

Version  1.0 


Kevin  B.  Smith 
Marine  Physical  Laboratory 
Scripps  Institution  of  Oceanography 
University  of  California,  San  Diego 

Frederick  D.  Tappert 
Applied  Marine  Physics 

Rosenstiel  School  of  Marine  and  Atmospheric  Science 
University  of  Miami 


MPL  Technical  Memorandum  432 


DT1C 

electe 
OCT  141993 

E 


D 


May  1993 


Approved  for  public  rdeue:  distribution  unlimited 


University  of  California,  San  Diego 
Scripps  Institution  of  Oceanography 


93  i0  6  1 5  5 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporing  burden  (or  this  collection  of  information  is  estimated  to  average  i  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  date  sources 
gathering  and  maintaining  the  data  need  ed.  and  com  pie  ting  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of 
this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway.  Suite  1204  Arlington,  VA  22202-4302.  and  to  the  Office  of  Management  and  Budget  Paperwork  Reduction  Project  (0704-0188).  Washington  DC  20503 


1.  Agency  Use  Only  (Leave  Blank). 


4.  Title  and  Subtitle. 


3.  Report  Type  and  Dates  Covered. 

Technical  Memorandum 


UMPE;  The  University  of  Miami  Parabolic  Equation  Model  (Version  1.0) 


5.  Funding  Numbers. 

N0001 4-93-1 -0062 


6.  Author(s). 

Kevin  B.  Smith  and  Frederick  D.  Tappert 


7.  Performing  Monitoring  Agency  Names(s)  and  Addresses). 

University  of  California,  San  Diego 
Marine  Physical  Laboratory 
Scripps  Institution  of  Oceanography 
San  Diego,  California  92152 


Project  No. 
Task  No. 


Report  Number. 


MPL-U-27/93 


9.  Sponsoring/Monitoring  Agency  Name(s)  and  Addresses). 

Office  of  Naval  Research 
Department  of  the  Navy 
800  North  Quincy  Street 
Arlington,  VA  22217-5000 
Code  1 1 250 A 


10.  ^Kmonji^/lgonitoring  Agency 


12a.  Distribution/Availability  Statement. 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  Distribution  Code. 


13.  Abstract  (Maximum  200  words). 

Solving  the  acoustic  wave  equation  using  the  parabolic  approximation  is  a  popular  approach  of  many  available 
underwater  acoustic  models.  Here  we  develop  and  present  a  verion  of  the  PE  model  developed  at  the 
University  of  Miami  from  1989  to  the  present  under  the  guidance  of  Professor  Fred  Tappert. 


14.  Subject  Terms. 

acoustic  wave  equation,  parabolic  approximation,  volume  attenuation, 
broadband  travel  time  analysis 


15.  Number  of  Pages. 

93 


16.  Price  Code. 


19.  Security  Classification 
of  Abstract.. 

Unclassified 


20.  Limitation  of  Abstract. 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Std.  239-18  298-102 


Contents 


1.  General  theory  1 

2.  Boundary  conditions  1 S 

2.1  Water/bottom  interface  16 

2.1.1  Sound  speed  discontinuity  16 

2.1.2  Density  discontinuity  19 

2.2  Surface  interface  25 

3.  Volume  attenuation  26 

3.1  Water  volume  attenuation  27 

3.2  Bottom  volume  attenuation  28 

3.3  Effective  attenuation  due  to  shear  28 

4.  Rough  interface  forward  scatter  30 

4.1  Water/bottom  interface  33 

4.2  Surface  interface  33 

4.2.1  Exact  forward  scatter  34 

4.2.2  Approximate  forward  scatter  37 

4.3  Effect  of  near-surface  bubbles  39 

5.  Broadband  travel  time  analysis  41 

6.  Acoustic  particle  velocity  44 

7.  Source  functions  47 

8.  Filters  and  sponges  54 

9.  Calculations  of  transmission  loss  55 

10.  Future  upgrades  57 

10.1  Automatic  cq  selection  58 

10.2  Azimuthal  coupling  59 

11.  Numerical  implementation  and  organization  63 

12.  Examples  71 


WIC  quality  INSPECTED  2 


UMPE: 

The  University  of  Miami  Parabolic 

Version  1.0 


Kevin  B.  Smith 
Marine  Physical  Laboratory 
Scripps  Institution  of  Oceanography 
University  of  California,  San  Diego 

Frederick  D.  Tappert 
Applied  Marine  Physics 

Rosenstiel  School  of  Marine  and  Atmospheric  Science 
University  of  Miami 

MPL  Technical  Memorandum  432 
May  1993 


Approved  for  public  release;  distribution  unlimited 


Contents 

1.  General  theory  1 

2.  Boundary  conditions  15 

2.1  Water/bottom  interface  16 

2.1.1  Sound  speed  discontinuity  1 6 

2.1.2  Density  discontinuity  19 

2.2  Surface  interface  25 

3.  Volume  attenuation  26 

3.1  Water  volume  attenuation  27 

3.2  Bottom  volume  attenuation  28 

3.3  Effective  attenuation  due  to  shear  28 

4.  Rough  interface  forward  scatter  30 

4.1  Water/bottom  interface  33 

4.2  Surface  interface  33 

4.2.1  Exact  forward  scatter  34 

4.2.2  Approximate  forward  scatter  37 

4.3  Effect  of  near-surface  bubbles  39 

5.  Broadband  travel  time  analysis  41 

6.  Acoustic  particle  velocity  44 

7.  Source  functions  47 

8.  Filters  and  sponges  54 

9.  Calculations  of  transmission  loss  55 

10.  Future  upgrades  57 

10.1  Automatic  c0  selection  58 

10.2  Azimuthal  coupling  59 

1 1 .  Numerical  implementation  and  organization  63 

12.  Examples  71 


Solving  the  acoustic  wave  equation  using  the  parabolic  approximation  is  a  popular 
approach  of  many  available  underwater  acoustic  models.  Here  we  develop  and  present  a  version 
of  the  PE  model  developed  at  the  University  of  Miami  from  1 989  to  the  present  under  the  guidance 
of  Professor  Fred  Tappert.  The  model  has  aptly  been  named  the  UMPE  model.  Fundamentally  a 
research  model,  the  numerical  approaches  used  here  may  be  compared  to  other  PE  models  and, 
subsequently,  may  begin  a  discussion  on  the  diversity  and  validity  of  current  PE  methods. 

This  report  is  based  on  a  collection  of  published  articles  by  various  authors  and  a  compila¬ 
tion  of  unpublished  lecture  notes  given  by  Fred  Tappert.  The  motivation  here  is  an  attempt  to 
address  the  fundamental  properties  of  PE  modeling  and  the  implementation  within  the  UMPE 
model.  The  framework  for  the  algorithms  and  a  description  of  the  source  code  implementation  will 
also  be  given.  As  upgrades  to  both  the  code  and  this  text  will  occur  inevitably,  we  are  defining  this 
model  as  version  1.0.  Users  are  encouraged  to  contact  either  author  with  questions  or  numerical 
problems  (i.e.,  bugs). 


1.  General  theory  of  PE  approximations 


We  begin  by  representing  the  time-harmonic  acoustic  field  in  a  cylindrical  coordinate  sys¬ 
tem  by 

P(r,z,  <p,  id/)  -  p(r,z,  q>)  e~l(al  (1.1) 

Substituting  this  into  the  wave  equation  in  cylindrical  coordinates  leads  to  the  Helmholtz  equation, 


+  ~2^r  +  ^4  +  kW  (r’  *  «P )p  -  -4KP0i>  (x  ~  xs) 

r  dr  dr  r *  dz" 


0)  co 

where  kQ  -  —  is  the  reference  wavenumber,  n(r,  z,  (p)  -  — - — r-  is  the  acoustic  index  of 

c o  c(r,  z,  cp) 

refraction,  c0  is  the  reference  sound  speed,  and  c  (r,  z,  cp)  is  the  acoustic  sound  speed.  It  is  within 

c(r,z,  cp)  that  all  features  of  the  environment  are  represented.  The  source  function  is  that  of  a 
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point  source  at  coordinates  (r  -  0,  z  -  zs)  with  reference  source  level  PQ  defined  as  the  pressure 
amplitude  at  a  reference  distance  of  R0  -  1  m,  and 

6(i)  -  ^6(z-rs)  6(r)  (1.3) 

Primarily,  acoustic  energy  propagates  outward  from  acoustic  sources  in  the  ocean  in  the 
horizontal  direction.  We  therefore  represent  the  pressure  field  by  an  outgoing  Hankel  function 
slowly  modulated  by  an  envelope  function, 

/>(r,r,<p)  -  ^(r.z.rp)/^0  (k0r)  .  (1.4) 

In  the  far-field,  the  Hankel  function  can  be  approximated  by 

•  (15) 

thus  an  alternative  relationship  between  the  acoustic  pressure,  p  (r,  z,  <p) ,  and  die  slowly  modulat¬ 
ing  envelope  function,  ip  ( r ,  z,  cp) ,  is 


P  ( r ,  s,  q>)  -  P0  JyQ  (r,  z ,  q>)  e'*°r  .  (1 .6) 

This  is  the  standard  definition  of  the  so-called  “PE  field  function”  ip  scaled  such  that  at  r  -  R0, 
|ip|  -  1  and  \p\  -  PQ.  Substituting  (1.6)  into  (1.2)  and  dropping  the  source  contribution  to  the 
far-field  solution  yields  the  defining  equation  for  the  evolution  of  the  PE  field  function. 


J?  +  2,k°Tr  + 


1  jjp 

r2d( p2 


,2/2  , 

— r-  +  kn  In  -  1  + 

a*-2  \ 


4  i&r) 


^  Ip  - 


(1.7) 


Because  this  is  a  far-field  approximation,  we  choose  to  neglect  the  two  terms  containing  the  4r  fac¬ 


tors.  This  will  automatically  invoke  the  uncoupled  azimuth  (UNCA)  approximation  by  dropping 

the  — j  term.  Note  that  this  is  strictly  valid  only  for  — « ip .  Should  the  environmental  condi- 
d<p  3<p“ 
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•» 

d“il> 

tions  become  such  that  — -  »  ip  and  is  non-negiigible  with  respect  to  then  azimuthal 

dq>“  °r 

coupling  should  not  be  ignored.  We  will  return  to  this  topic  near  the  end  of  this  report. 

At  this  point,  the  only  approximations  we  have  applied  to  the  wave  equation  are  a  far-field 

assumption  and  little  or  no  azimuthal  variations.  Implicit  in  the  far-held  analysis  is  the  small  angle 

approximation,  evident  from  the  expression  for  the  uniform  ocean  Green's  function 


1 


1/2 


i  a 


[r2+(z-zs)2] 


where  zs  is  the  source  depth.  For  small  angles  of  propagation,  |6|  «* 


(1.8) 


«  1,  so 


p»-e 

r 


1 ... /V 


-  -Fife 
Jr 


(1.9) 


The  corresponding  form  for  i|>  is  found  to  satisfy  a  parabolic  type  equation.  The  general  form  of 
the  parabolic  approximation  to  the  wave  equation  simply  follows  then  from  the  acknowledgment 

that  tp  is  a  slowly  varying  function  in  range  and  we  may  neglect  the  term  —j .  Thus,  Eq.  (1.7) 


dr 


takes  the  form 


dtp 


i  d  ik  o  ,2  n 
+  — (w  -  1)^ 


(1.10) 


dr  2kodJ 

Eq.  (1 . 10)  was  first  introduced  to  the  underwater  acoustics  community  by  Tappert  (1974). 

We  have  now  reduced  a  second  order  differential  equation  to  a  first  order  one,  thereby 
allowing  solutions  via  a  non-iterative  marching  algorithm.  Note  that  we  may  rewrite  Eq.  (1 . 10)  as 


lko'rr  -  -  (T«r*uop>  t 


(1.11) 


where  the  operators 


op 


_ L_£i 

dz  2*o  dz2 


(1.12) 
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and 


Uop  -  U(r,  :,tt)  -  -i  (n2-  1)  (1.13) 

This  representation  of  the  operators  as  kinetic  and  potential  energy  operators  is  especially  insight¬ 
ful  when  one  wishes  to  form  the  ray  equations  which  have  Hamiltonian  form  (e  g.,  Smith  et  al., 
1992).  In  Eq.  (1.11),  the  function  ip  is  a  vector  (in  z)  in  Hilbert  space.  The  relationship  between 
the  values  of  ip  at  different  ranges  can  now  be  expressed  as 

ip(r  +  Ar)  -  d>(r)ip(r)  (1.14) 

To  propagate  the  solution  out  in  range  requires  a  representation  of  the  propagator  <2>  (r) 

There  are  three  common  methods  of  computing  PE  solutions:  (1 )  the  split-step  Fourier  (PE/ 
SSF)  method  (Hardin  and  Tappert,  1973),  (2)  the  implicit  finite  difference  (IFD-PE)  method  (e  g., 
Lee  and  Botseas,  1982),  and  (3)  the  finite  element  (FEPE)  method  (e  g.,  Collins,  1988).  Since  the 
UMPE  model  uses  the  first  technique,  we  shall  isolate  our  discussion  to  the  implementation  of  the 
PE/SSF  method.  This  is  easily  accomplished  by  approximating  the  propagator  function  by 

4>(r)  ~e~ik°Bopir) Ar  (1.15) 

where 

H„p(r)  -  L^^dr'Hop(r')  .  (1.16) 

The  validity  of  this  approximation  can  be  seen  from  the  following  argument.  Suppose  we  divide 
Eq.  (1.11)  by  ip  to  obtain 

JL(/nip)  -  -ikQHop(r)  .  (1.17) 

Integrating  this  yields 

//np  (r  +  Ar)  -  Im p  (r)  -  ikQfr  *  ^ dr' Hop  (r') 
or 


ip(r  + Ar) 


-ikQJr*Srdr'Hop(r') 


ip(r) 


(1.18) 
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This  is  exactly  the  equation  we  developed  in  ( 1 . 1 4),  ( 1 . 1 5),  and  (1.16).  It  is  not  formally  valid, 
however,  since  Hop  is  an  operator  and  we  cannot  divide  (1.11)  by  tp.  In  other  words,  it  is  only  an 


approximation  because 


dHn„  -i 


L  dr 


The  formal  solution,  using  a  Dyson  time  evolution 


operator  (Sakurai,  1985),  would  be 

rr  dr'H^(r') 

i|>(r  +  Ar)  -  Te  Jr  t|>(r)  (1.19) 

where 


+  Ar 


dr'jrdr”Ht 


+  ( higher  order  terms )  .  ( 1 .20) 

Thus,  from  Eq  .  (1 .20),  we  can  evaluate  the  first  order  correction  to  Eq.  ( 1 . 1 8). 

Assuming  (1 . 1 8)  to  be  valid,  we  still  must  evaluate  Hop  defined  in  Eq.  (1 . 16).  Two  com¬ 
mon  approximations  are 


Ho,  -  H  lr+ltir) 

r 


0  21) 


and  simply 

Hop  -  Hop(r)  .  (1.22) 

These  are  sometimes  referred  to  as  the  “centered”  and  “end  point”  schemes,  respectively.  The 
interpretation  of  these  approximations  is  that,  over  the  range  step  r  tor  +  Ar,  the  operator  (hence 
the  environment)  is  sampled  at  either  the  middle  or  the  beginning  of  the  range  step.  Presumably, 
if  Ar  is  small  enough  the  differences  between  the  solutions  are  negligible.  The  UMPE  model 
implements  Eq.  (1.22)  and  we  shall  isolate  the  remainder  of  our  discussion  accordingly. 

The  operator  Uop  is  simply  a  multiplication  operator  in  z-space  and,  hence,  is  a  diagonal 
matrix.  The  operator  Top  is  not  diagonal  in  z-space  so  different  depth  eigenfunctions  are  coupled. 

In  wavenumber  space,  however,  the  corresponding  operator  Top  is  diagonal.  It  is  desirable,  there- 
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fore,  to  separate  the  application  of  each  operator,  one  in  r-space  and  one  in  A-space  Using  the 
Baker-Campbell-Hausdorff  expansion  (Bellman,  1964),  we  may  write 


A+B  A  BJA,  B)  *  [A,  +  [B,  ♦. 

►  ®  V  v 


(123) 


where  A  -  -ik^rTop  and  B  -  -ik0ArUop  Since  both  Top  and  Uop  are  smp’i  then  we  assume 


their  products  are  of  second  order  and  negligible.  Finally  then,  we  have 

<p(r)  _  e-ik<A'T*e-ik'£rU., 


(1.24) 


Note  that  this  separation  of  Hop  into  two  components,  each  of  which  is  diagonal  in  some  represen¬ 
tation  and  can  be  applied  independently  of  the  other,  is  presumed  by  application  of  the  SSF 
integration  scheme.  The  various  approximations  used  to  separate  the  operator  Hop  are  typically 
used  to  distinguish  one  type  of  PE/SSF  model  from  another. 

Note  from  Eq.  (1.24)  that  if  there  are  no  losses  present  (i.e.  ImT  -  ImU  -  0)  then 


IIOMII  -  1  • 

and  <J>  (r)  is  a  unitary  operator.  Therefore,  the  normalization  condition  is 


(125) 


H(r)\\  -  f\y(r,z)\2dz  -  constant 


(1.26) 


In  other  words,  because  of  the  formulation  of  the  propagator,  the  PE/SSF  scheme  is  conservative. 
There  are  no  intrinsic  losses  due  to  the  numerical  scheme. 

The  general  algorithm  behind  the  PE/SSF  implementation  is  ^hen  as  follows.  The  PE  field 
function  r|>  is  specified  at  some  range  r  in  the  r-domain.  A  transformation  is  made  to  the  ^-domain 

followed  by  a  multiplication  of  the  A-space  rator  e~,k^rTop  The  result  is  then  transformed  again 

to  the  2-domain  and  is  followed  by  a  multiplication  of  the  r-space  operator  e~'*oArtV  (The  actual 
order  is  irrelevant  since  the  commutation  of  these  operators  is  considered  insignificant  but  this  is 
the  order  imposed  in  UMPE.)  The  final  result  is  the  field  function  at  r  +  Ar.  The  FFT  subroutine 
employed  in  the  numerical  code  assumes  the  convention 
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V(-)  -  FFT  (ty  {k) ) 


(1.27) 


♦  («  -  [FFT  (ip*  (--))]*  . 
Therefore,  the  PE/SSF  implementation  can  be  represented  by 


(1.28) 


V  (r  +  At,  r)  -  x  FfT  ^-«*0Arfv(r,*)  x  [pFT  (/%  r)  j  ]  *  j  f  (i.29) 


where,  in  A-space, 


**«<*)  -i<f)  • 


(1  30) 


Previously,  we  have  assumed  the  operators  took  the  forms  defined  by  Eqs.  (1.12)  and 
(1.13).  These  forms,  which  followed  from  the  derivation  of  the  parabolic  equation  (1 . 10),  are  com¬ 
monly  referred  to  as  the  “standard  PE”  or  SPE  forms  and  are  only  one  set  of  a  number  of  various 
operator  forms.  To  obtain  other,  higher  order  forms,  we  return  to  the  original  wave  equation  (1.2). 
Still  ignoring  the  source  term  and  the  azimuthal  coupling  term,  we  now  define  the  pressure  field  as 


p(r,z)  -  -=u(r,z) 

Jr 


(1.31) 


The  function  u  ( r ,  z)  is  identical  to  the  pressure  field  in  two  dimensions  and  the  term  -=.  accounts 

Jr 

for  azimuthal  spreading.  Substituting  (1.31)  into  the  Helmholtz  equation  in  two  dimensions  yields 


the  far-field  UNCA  expression 


2  ^ 

d  U  d~u  , 2  2,  x  A 
—  +— T+kon  \r,z)u  -  0 

dr-  dz2 


(1.32) 


We  introduce  the  operators 


P  -  — 
°p  dr 


(1.33) 


(1.34) 
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Eq.  (1 .32)  then  becomes 


-  o 


(1.35) 


which  can  be  factored  as 

(Pop  +  ikoQop)  ( Pop-ikoQop )  U  +  ik0  lPop>  Qop ]  u  -  0  •  o  36) 

The  commutator  [P0^  Q  }  is  assumed  negligible  and  is,  in  fact,  exactly  zero  in  layered  media. 
Eq.  (1.36)  therefore  represents  the  combination  of  incoming  and  outgoing  waves.  The  outgoing 
wave  satisfies 

Popu  “  iko QopU  (137) 

or 

-  QoP“  0  38) 

Eq.  (1.38)  is  observed  to  have  a  form  similar  to  Eq.  (1.11).  The  trick  is  now  to  develop  an  approx¬ 
imation  for  Qop  that  separates  into  a  r-space  operator  and  a  £-space  operator. 

The  UMPE  model  allows  the  user  to  choose  from  five  different  Qop  approximations  which 
we  shall  now  derive.  We  begin  by  introducing  the  notation 

e  -  n2  -  1  (1.39) 

and 

1  a2 

H  0- *0) 

*5  a--2 

so 

Qop  -  (n  +  t  +  Dw  (141) 

The  first  approximation  follows  from  the  assumption  that  both  e  and  p  are  small  compared  to  unity. 
A  binomial  expansion  then  yields 

Qop~Qi  -  5K+5S+I  (1  42) 
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Substitution  of  (1.39)  and  (1.40)  shows  that  this  is,  within  an  additive  constant,  the  standard  PE 
operator. 


2,  -  Qspe  -  ^-4*  ^ <"2  - 1)  + 1  (1.43) 

2*o  dz  1 

ic  2 

Interpreting  the  operator  p  as  —  -  sin  0  indicates  that  the  assumption  p  «  1  implies 

2  2 

sin  0  •»  0  «  1 ,  hence  this  is  a  small  angle  approximation.  The  condition  e  «  1  is  simply  inter¬ 
preted  as  assuming  primarily  stratified  media  as  is  typical  of  most  ocean  regions.  In  fact,  since 

e  -  ( n  —  1 )  «  1 ,  we  can  further  approximate  this  by  e  -  2  («  -  1)  to  obtain 

Q,  -  -L-^+(»-l)  +  l  .  (1.44) 

2  kid;2 

Both  of  these  approximations  were  first  recognized  by  Tappert  (1977)  in  his  original  Springer- Ver- 
lag  article. 

To  obtain  the  customary  formula  for  the  PE  field  function  ip,  we  first  give  it  the  usual  enve¬ 
lope  definition 

u(r,z)  -  )(r,z)ek°r  .  (1.45) 

Substituting  this  into  (1.38)  yields 

—  -  -  ik0q  +  ik0Qopy  .  (1.46) 

Application  of  Qop  -  Qspe  produces  exactly  Eq.  (1.10).  Comparing  (1.46)  with  (1.11)  indicates 
the  simple  relationship 

Qop "  1  -Hop  -  1  -(Top+Uop)  (1.47) 

which  easily  allows  a  split-step  Fourier  implementation.  Equivalently  then,  our  goal  is  to  develop 
expressions  for  Top  and  Uop  based  on  approximations  of  Qop. 

A  better  approximation  to  Qop,  also  introduced  by  Tappert  (1977),  assumes  that  only  the 
sound  speed  variations  are  small,  i.e.  e  «  1 .  We  may  then  write 
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(148) 


Qi  -  (1  +n)1/!  +  ^e 

The  restriction  on  propagation  angles  has  now  been  removed.  In  fact,  for  homogeneous  media 
(e  -  1 ) ,  Eq.  (1 .48)  is  observed  to  be  an  exact  representation  of  the  original  Helmholtz  operator. 

A  higher  order  approximation  introduced  by  Thomson  and  Chapman  (1983)  is  actually  a 
combination  of  Eqs.  (1 .44)  and  (1 .48).  It  is  based  on  an  operator  splitting  by  Feit  and  Fleck  (1978) 
which  is  formally  valid  only  when  e  and  p  commute.  Commonly  referred  to  as  the  “wide-angle” 
approximation  (WAPE),  it  has  the  form 

Qa-Qwape-  0*|*),/1+[<1+«)1/I-1]  (149) 

Invoking  the  operator  identity 

(1+p)1^  -  l+iitO  +  lO^+l]"1  (1.50) 

and  formulating  Eq.  (1 .49)  in  terms  of  Top  and  Uop  as  in  Eq.  (1 .47)  leads  to 

2f  2 1/1  r1 

W  -  -TJ  O  +  1  0  51) 

dz  L  dz  J 

and 

<W  -  •  0.52) 

In  wavenumber  space,  we  may  express  (1.51)  as 

,r  2  1/2  I-1  r 

TwapeW  -  (£>'  i  +1  -  1  -  1  -  <r>  053) 

*o  L  \  *o  /  J  L  . 

Note  that  modes  with  k  >  k0  are  evanescent  since 

r  .  2  -i1'2 

WfHJ-i-Mf)-1  (153’) 

.  *0  J 

The  final  approximation  we  shall  consider  was  introduced  by  Berman  et  al.  (1989)  and  is 
referred  to  as  the  “modified  LOGPE”  or  simply  LOGPE.  It  follows  from  assuming  the  potential 
energy  function  Uop  has  the  form 
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U LOG  PE  “  ~ln(n)  0  54) 

Comparing  the  ray  equations  derived  from  the  Helmholtz  equation  with  the  general  formulation  of 
the  ray  equpt'ons  in  terms  of  Top  and  Uop  leads  to  the  definition  of  the  kinetic  eneigy  term  consistent 
with  (!.  4) 


T logpe  ~  -/»£cos(-j^) 
In  the  wavenumber  domain,  this  becomes  simply 


(1.55) 


(1.55’) 


We  note  that  for  small  propagation  angles  ( A  «  1 ) , 

ko 


LOGPE 


j! 


2 


-  T, 


SPE 


while  for  nearly  uniform  media  ( n  «*  1 ) , 


(1.56) 


^LOGPEm~2  (n2  ~  ^  "  ^SPE  (157) 

Thus,  LOGPE  reduces  to  SPE  in  situations  where  the  latter  approximation  is  valid.  Unfortunately, 
*  k  3t 

the  operator  TLOgpe( *)  's  undefined  for  j-  je  - .  To  avoid  this,  a  sine-squared  taper  function  is 

*0  2 


1  In  3i  k  3i 

applied  over  the  outer  -  range  of  —  to  - .  For  7-  >  -,  this  function  is  set  to  zero. 

O  10*.  Ke\  2 


Of  the  five  PE  approximations  described  so  far,  the  last  two  are  expected  to  yield  the  most 
accurate  results.  The  most  common  PE  implementation  currently  is  the  WAPE,  and  this  is  the  ver¬ 
sion  used  in  the  Navy  standard  model  (Holmes  and  Gainey,  1991).  However,  the  WAPE  and  other 
so-called  “higher-order”  models  are  still  not  exact  and  may  occasionally  produce  results  that  are 
worse  than  predicted  by  the  SPE.  The  most  famous  example  of  this,  sometimes  referred  to  as  the 
Porter  duct  problem,  was  defined  in  Test  Case  7  of  the  PE  Workshop  II  (Chin-Bing  et  al.,  1993). 
In  such  instances,  it  is  usually  found  that  the  error  results  from  extra-sensitivity  to  the  choice  of 
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reference  sound  speed.  In  fact,  the  choice  of  c0  is  the  one  ambiguous  feature  of  all  PE  models,  and 
till  now  we  have  ignored  this  effect.  Later  in  this  report,  a  scheme  for  computing  a  default  value 
for  c0  will  be  discussed.  However,  modelers  are  hereby  forewarned  that  no  method  is  foolproof. 
The  best  approach  is  to  vary  c0  and  look  for  fluctuations  in  the  calculations  Currently,  the  UMPE 
model  only  takes  user  input  cq  values. 

Because  of  the  ever  present  ambiguity  in  the  selection  of  cq,  it  is  highly  desirable  to  develop 
a  model  that  is  c0-insensitive  such  that  significant  changes  in  the  choice  of  c0  will  not  affect  the 
final  result.  Tappert  (1991b)  developed  a  rigorous  definition  of  such  a  PE  model.  Application  of 
the  resulting  code  was  found  to  eliminate  the  sensitivity  in  the  Porter  duct  problem  However, 
implementation  was  complicated  and  required  a  transformation  of  the  function  ip  (z)  to  a  new 
function  ip  (z)  in  a  transformed  space.  In  a  range-dependent  environment,  such  a  transformation 
would  be  required  at  each  range  step  thereby  greatly  increasing  the  run  time.  Since  one  of  the  main 
advantages  of  the  PE/SSF  code  is  the  speed  with  which  the  acoustic  field  can  be  computed,  a  cQ- 
insensitive  version  has  not  been  implemented  in  the  UMPE  code.  Again,  the  user  is  reminded  to 
beware  of  these  highly  sensitive  (but  uncommon)  problems. 

Two  other  ambiguous  variables  must  also  be  introduced  in  the  numerical  implementation. 
As  in  all  models,  a  discretization  of  the  environment  is  required  and  defined  by  the  mesh  size 
(A r,  Az) .  Because  the  depth  mesh  influences  the  wavenumber  increments  M,  we  may  define  a 
default  value  for  Az,  hence  the  transform  size  N,  by  considering  a  lower  limit  on  allowable  angles 
of  propagation.  Since  jV  wavenumber  values  will  be  spread  over  the  range  +kmaxlo  -kntaxy  it  follows 
that 


.  N\t  At  23t 

*max  m  2  ”  -j. 


where  Zj  is  the  total  computational  depth,  so 


(1  58) 


Furthermore,  the  wavenumbers  are  related  to  the  angles  of  propagation  by 

k  -  £0sin6  . 


(1.60) 


It  follows  that  for  a  given  maximum  angle  of  propagation,  the  minimum  transform  size  required 


must  satisfy 


,v  ,^£ls,„0 


(1.61) 


To  define  an  upper  bound  on  the  range  step  size,  A rmax,  we  consider  the  analogy  of  phys¬ 
ical  optics.  In  the  vicinity  of  a  focus,  the  signal  will  vary  significantly  over  a  horizontal  range  of 


(1.62) 


where  f  is  the  f-number  of  the  focusing  lens. 


(1.63) 


R  is  the  focal  length  and  2B  is  the  effective  aperture.  In  the  underwater  acoustics  problem,  the  focal 
length  can  be  represented  by  the  distance  between  convergence  zones,  CZ,  and  the  effective  aper¬ 
ture  is  roughly  the  depth  of  the  ucean  or  approximately  ^zT  (the  factor  of  ^  for  the  depth  will  be 


explained  later).  Then  (1.63)  becomes 


CZ  1 


064) 


where  6  is  the  angle  of  propagation.  Combining  ( 1 .62)  and  ( 1 .64)  yields  the  upper  limit  on  foe 
range  step  size 


.s-  •  0  65) 

k0*m  *mox 

A  similar  analysis  suggests  that  the  maximum  vertical  mtsn  size  is  given  by 

2f  7 

Ar  y — 7^ —  (1 .66) 

""  *0  *0 

which  can  be  shown  to  yield  roughly  the  same  order  for  foe  transform  size  as  Eq.  (1 .61). 
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From  the  above  analysis,  it  is  obvious  that  if  a  particular  problem  is  known  to  contain  only 
small  angle  propagation,  the  mesh  size  ( A r,  A :)  may  be  increased  and,  subsequently,  the  run-time 
will  be  reduced.  Conversely,  for  problems  where  very  large  angle  propagation  is  expected  to  be 
important  a  small  mesh  size  may  be  required.  If  the  user  is  unsure,  zero  values  may  be  input  for 
the  depth  transform  size  and  the  range  step  size  which  will  prompt  the  model  to  use  the  default  val¬ 
ues  defined  by  Eqs.  (1 .6 1 )  and  (1 .65).  A  value  for  Bmax  is  then  defined  by  the  input  source  angular 

width  or  30°,  whichever  is  greater. 


We  will  now  review  the  main  results  of  this  section.  The  UMPE  model  can  solve  five  dif¬ 
ferent  PE  type  approximations  to  the  wave  equation.  The  general  form  of  each  equation  is 

dii) 

Yr  "  -'*o (ToP+Uoo)V 


'op' 


and  the  acoustic  pressure  is  then  defined  by 


0  ,  ik0r 


P  - 

The  numerical  algorithm  employed  is  the  split-step  Fourier,  or  SSF,  such  that  the  solution  for  ip  is 
marched  out  in  range  according  to 


tp  (r,  r  +  Ar)  -  fiko^uoM  r)  x  FpT  {  r)  x  [FFT  (ip*  (r,r))]*}  , 

where,  by  definition, 

*(*>  -  FFT(tj>(*)) 
and 


♦  W  -  [FFT (ip*  (r))]* 

The  five  types  of  PE  approximations  and  their  corresponding  operators  are 
1)  SPE: 

1  *  2 

TsPE(k)  -  i(^)  , 
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2)  Variant  SPE: 


3)  Variant  WAPE: 


4)  WAPE: 


5)  LOGPE: 


Tl'SPE(k)  “ 


Uyspsi2)  m  ~[w(^)  ~  1]  > 


TvwapeW  ■  1  “ 


i-(r> 

*0  J 

1  r„2, 


,1/2 


Uvwape(z)  ™  “j*  tw  (r)  ~  n  * 


l-(£) 

*0  J 


,1/2 


TwapeW  •  1  - 

UwAPE^*}  m  “  I]  * 


TlogpeW  m  -/w^cos  (^)  J  ’ 

^ LOGPE  (“)  “  —  /«  [  w  ( 2^)  ] 


2.  Boundary  conditions 

The  only  parameter  introduced  in  the  last  section  for  any  of  the  PE  approximations  was  the 
reference  sound  speed  c0.  In  a  sense,  this  parameter  is  related  to  the  accuracy  of  the  “small  angle” 

k 

approximation  since  it  will  serve  to  scale  the  relative  wavenumber  values  .  The  best  choice  of 


c0  should  then  produce  the  most  accurate  results  for  a  given  problem.  However,  in  the  numerical 


implementation  of  the  PE,  additional  parameters  must  be  introduced  to  define  the  discretely  sam¬ 
pled  environment  in  a  way  that  also  yields  the  best  results.  This  is  primarily  a  concern  in  regions 
where  boundaries  exist.  The  representation  of  these  boundaries  in  a  PE/SSF  implementation  is  the 
subject  of  this  section. 


2.1  Water/bottom  interface 


The  UMPE  model  treats  the  bottom  as  a  fluid  of  contrasting  sound  speed  and  density  from 
that  of  water.  In  addition,  the  UMPE  model  allows  for  an  additional  bottom  layer  to  exist  on  top 
of  the  basement  to  allow  for  effects  of  sediment  layers  to  be  included.  Within  either  bottom  vol¬ 
ume,  the  PE  environmental  potential  function,  Uop  (z) ,  is  defined  as  before  in  terms  of  the  local 


acoustic  index  of  refraction,  n  (r)  -  -j-r ,  where  c(z)  now  includes  the  sound  speed  within  the 

c{z) 

bottom  volume.  The  effect  of  approximating  the  bottom  as  a  fluid  is  the  neglect  of  shear  wave 
propagation.  When  the  true  bottom  does  support  shear  waves,  the  conversion  of  compressions! 
energy  incident  on  the  interface  into  downward  propagating  shear  energy  is  treated  as  a  loss.  In 
this  manner,  the  bottom  properties  are  replaced  by  equivalent  fluid  properties  that  produce  the  cor¬ 
rect  reflection  from  the  interface.  The  treatment  of  shear  as  a  loss  mechanism  is  discussed  further 
in  a  following  section. 


2.1.1  Sound  speed  discontinuity 

We  assume  the  interface  between  the  bottom  of  the  water  column  and  the  top  of  the  base¬ 
ment,  or  sediment,  layer  is  characterized  by  a  sharp  contrast  in  sound  speed.  In  a  numerical  code 
with  finite  sampling  and  recurrent  use  of  FFT’s,  it  is  desirable  to  use  smoothly  varying,  continuous 
functions  to  avoid  artificial  reflections,  aliasing,  and  noise  from  entering  into  the  calculation. 
Therefore,  we  seek  to  find  a  smooth,  continuous  function  of  variable  scale  which  can  accurately 
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reproduce  the  physical  effects  of  a  discontinuous  jump  in  sound  speed  at  the  water/bottom  inter* 
face. 


We  write  the  interface  condition  for  the  sound  speed  as 

c(i)  -  C,  (s)  +  [cb  (z)  -  c,  (r)  ]  H  (z  -  zb)  (2. 1) 

where  we  shall  assume  that  the  sound  speed  above  the  interface  at  z  -  zb  has  a  constant  value  of 
cw  and  below  the  interface  has  a  value  cb.  The  Heaviside  step  function  is  defined  by 

0,  £<0 


HiX>) 


£-0 


(2.2) 


1,  £>o 

where  t,mz-zb.  From  the  theory  of  generalized  functions  (Lighthill,  1958),  we  may  replace  H(t,) 
by  any  smooth  function  within  a  class  of  generalized  functions  that  produces  the  same  overall 
effect  (i.e.,  has  similar  moments). 

One  such  function  satisfying  the  above  criteria  involves  the  hyperbolic  tangent  function. 


H(K.) 


(23) 


or,  equivalently. 

Hi®  -  ( 1  +  e~**/L')  1  .  (2.3’) 

This  function  has  the  properties 

//(£)-  0,  £ «  0;  (2.4a) 

»(5)  -  i,  S  -  0;  (2.4b) 

and 

Hit,)  ->  1,  C  »  0  •  (2.4c) 

Furthermore,  the  derivative  of  the  Heaviside  function  is 

H'(K)  -  &(£)  (2.5) 

where  6(£)  is  the  Dirac-delta  function  and  is  characterized  by 
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1 


(26) 


pax 

-00 

Similarly,  the  derivative  o f  //(C)  is 

??(£>  -  6(0  -  ^j-sech2(^j-)  (2.7) 

and  it  is  easy  to  show  that 

00 

pax-i  (2.8) 

-00 

is  also  satisfied.  This  mixing  function  //(£)  is  parameterized  by  a  characteristic  mixing  length, 
Lc.  It  is  obvious  from  the  above  analysis  that 

lim  Hi®  -  H&)  (2.9) 

and 

lim  5(0  -  6(0  (2.9*) 

Lc  —  0 

This  limiting  equality  can  be  shown  to  hold  for  higher  derivatives  as  well. 

We  have  now  introduced  an  additional  parameter  into  our  model,  the  sound  speed  mixing 
length,  Lc.  This  can  be  adjusted  by  the  user  in  an  attempt  to  create  the  most  realistic  interface  con¬ 
dition  for  reflections  from  a  sound  speed  discontinuity.  The  UMPE  code  then  employs  Eq.  (2.3’) 
to  mix  the  sound  spiced  profiles  above  and  below  the  interface  (or  interfaces)  atr  -  zb.  Experience 
has  shown  that  the  most  accurate  results  are  gained  by  defining  Lc  as  a  fraction  of  an  acoustic  wave¬ 
length.  Experience  has  also  shown  that  the  sampling  in  depth  A z  need  not  necessarily  be  less  than 
a  wavelength  in  which  case  Lc  ~  X/5  would  result  in  a  sharp  jump  from  cw  to  cb  over  one  mesh  point. 
Ironically,  this  is  what  we  intended  to  avoid.  The  interpretation  of  this  apparent  contradiction  is 
that  a  discontinuous  jump  in  sound  speed,  hence  a  jump  in  the  PE  potential  function  U0f/z),  does 
not  introduce  much  error  and,  in  fact,  yields  the  best  representation  of  the  reflection  condition.  A 
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default  minimum  value  for  Lc  has  therefore  been  set  at  Z,Cmm  -  —  Alternatively,  and  perhaps 

more  accurately,  the  index  of  refraction  could  be  smoothed  over  the  interface  Experience  has 
again  shown  that  this  produces  little  difference  in  the  final  result. 


2.1.2  Density  discontinuity 

The  effect  of  density  on  acoustic  propagation  has  not  yet  been  considered.  In  fact,  the  vari¬ 
ation  of  density  was  ignored  in  the  original  form  of  the  Helmholtz  equation.  In  a  fluid  with  a 
variable  density  p,  the  correct  form  of  the  Helmholtz  equation  is 

pV.(Ivp)+*5n2p  -  0  .  (2.10) 

We  make  the  substitution 


to  obtain  the  equation  defining  the  propagation  of  q. 


(2.11) 


(2.12) 


where  »'  is  the  “effective”  index  of  refraction  given  by 


n'2  -  n  +  — ^ 

2*o  ^ 


iv2P-|(ivp)' 


(2.13) 


Therefore,  we  may  solve  for  the  pressure  field  p(r,z),  now  defined  by 


(  \  D  P^° .  /  \  *V 

p{r,:)  -  P0  —y(r,z)e  , 

N  Por 


(2.14) 


by  marching  the  solution  of  the  PE  function  ty(r,r)  out  in  range  with  the  definition 

Uop(z)  -  Ux(z)+U2(z)  (2.15) 

where  f/j(r)  is  the  same  environmental  potential  function  previously  defined  and  U2(z)  accounts 
for  the  effect  of  the  density  discontinuity.  For  example,  if  the  SPE  approximation  is  employed,  one 
can  easily  show  that 
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where  we  have  assumed  p  -  p  (z)  only. 


The  reference  density  is  defined  as  that  of  pure  water,  p0  -  1  g/cm 3 ,  and  inputs  to  the 


p 

model  for  density  are  given  as  ratios,  i.e.  — .  The  UMPE  model  assumes  each  volume  (water,  sed- 

Po 

iment,  basement)  is  characterized  by  a  constant  density  profile.  For  example,  with  a  single  water/ 
bottom  interface, 

P(*)  -  Pw+ (P*-P„) #(---*)  (2  17) 

where  //(£)  is  the  Heaviside  step  function  described  previously.  Obviously,  the  function  U2(z)  is 
non-zero  only  in  the  vicinity  of  the  interface.  As  before,  we  wish  to  spread  this  discontinuity  over 
some  finite  region  in  terms  of  smooth  generalized  functions.  This  is  a  more  critical  problem  than 
before  because  U2(z)  depends  not  on  the  density  but  on  the  derivatives  of  the  density.  In  most  cur¬ 
rent  versions  of  PE/SSF,  use  is  made  of  Eq.  (2.16)  to  define  the  density  potential  function  and  the 
hyperbolic  mixing  function  Eq.  (2.3'),  with  its  subsequent  derivatives,  to  define  the  interface  con¬ 
dition.  However,  the  choice  of  the  density  mixing  length  Lp  has  always  been  somewhat  ambiguous 
in  this  definition.  If  Lp  is  chosen  too  laige  or  too  small,  regardless  of  how  well  the  function  is  sam¬ 
pled  in  space,  the  result  will  not  be  correct.  Recently,  Tappert  (1991a)  has  suggested  that  a  better 
form  for  this  function  is 


U2(z) 


where 


(218) 


'1-(PA)1/2' 

j  +  (pypb)'/2_ 


(2.19) 


For  small  density  contrasts,  this  is  equivalent  to  neglecting  the  second  term  in  Eq.  (2. 1 6).  The  main 
argument  used  to  justify  this  approximation  is  that  because  we  are  using  generalized  functions  to 


represent  the  density  discontinuity,  the  function  (A(r)  must  also  be  defined  in  terms  of  generalized 
functions  However,  the  second  term  in  (2. 16)  contains  the  square  of  a  generalized  function,  spe- 

_2 

cifically  6  {:')  .which  is  not  a  generalized  function.  A  detailed  analysis  by  Tappert(  199  Id)  which 
attempted  to  remove  as  much  singularity  as  possible  from  the  solution  in  the  vicinity  of  the  density 
discontinuity  showed  that  Eq.  (2. 18)  is  the  best  expression  to  use  in  the  SPE  approximation.  Fur¬ 
thermore,  this  formulation  inherently  produces  the  best  results  when  Zp  is  minimized,  i.e.  as 

Zp  -*  0  Note  that  Lp  must  still  be  large  enough  such  that  the  finite  depth  mesh  adequately  samples 
the  function  ^(r)  (which  is  not  a  simple  jump  as  U\  was).  We  shall  assume  (2  .18)  can  be  applied 
for  any  of  the  PE  approximations  used.  The  main  justification  for  this  is  the  expectation  that  most 
bottom  interacting  energy  will  be  near  grazing,  hence  small  angle  reflection  dominates  at  long 
ranges. 

We  now  examine  two  forms  of  mixing  function  available  in  the  UMPE  model  to  compute 
the  density  potential  function  fA(c): 

a)  The  first  approximation  of  the  Heaviside  step  function  is  the  same  hyperbolic  tangent 
function  used  to  mix  the  sound  speeds  across  the  interface,  i.e. 

H(t,)  -  (l+e'^V1  •  (2.20) 

The  first  derivative  is 

i 

H'iK,)  -  M£)  “  7- - - - *  (2.21) 

and  the  second  derivative  is 

</L  </L 

—  -  1  e  ( 1  —  e  ) 

H"iX,)  -  8'(C)  -  -4 - - - T2  ■  <2-22) 

Li  (  y  3 

p  (1+e  p) 

Substituting  this  into  the  density  potential  function  (2.18)  yields  one  of  the  formulas  used  in  the 
UMPE  model  to  compute  the  effect  of  a  density  discontinuity. 
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The  main  advantage  of  the  hyperbolic  smoothing  function  is  that  it  is  a  Cx  smooth  func 


tion,  i.e.  all  derivatives  are  continuous  everywhere.  The  main  disadvantage  of  this  type  of 
smoothing  function  is  that  it  has  “infinite”  extent.  Although  it  may  quickly  approach  a  negligible 
value  (even  zero  within  computer  precision)  away  from  the  interface,  it  may  still  extend  into  unde¬ 
sirable  regions. 


One  can  easily  show  that  the  extrema  of  6'  (£)  occurat£ex  -  Zpln(2±j/3)  «±1.32  Zp. 


At  these  values,  6'  (£  )  *-  ±0,  ~6.  Because  it  is  presumed  that  f/2(-)  yields  the  most  accurate 

L~ 

representation  for  small  Ip,  for  a  given  mesh  size  Ac  it  follows  that  the  optimum  choice  of  mixing 

At 

length  for  the  hyperbolic  tangent  function  is  approximately  .  This  produces  a  function  with 

extrema  exactly  one  depth  mesh  above  and  below  the  interface.  For  this  to  be  properly  sampled, 
however,  the  interface  must  also  lie  exactly  on  a  mesh  point.  To  accommodate  these  restrictions. 


a  minimum  value  of  Ln 

Pmin 


At 
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is  defined.  When  the  user  input  mixing  length  is  equal  to  or 


less  than  this,  the  potential  function  f/2(c)  is  centered  on  the  mesh  point  nearest  the  true  interface 
depth.  In  this  manner,  a  default  condition  is  set  if  the  user  inputs  Ip  =  0.  While  this  should  produce 
the  most  accurate  representation  of  U2(z),  it  may  introduce  an  error  in  the  depth  of  the  interface  by 


as  much  as  ^  Ac.  This  error  can  be  minimized  and  convergence  tested  by  decreasing  Ac  (i.e., 
increasing  the  transform  size). 

Experience  suggests  that  the  most  accurate  results  are  achieved  by  requiring  Zp  s  X.  Con¬ 
sequently,  this  implies  a  depth  mesh  limit  of  Armax  »  1 .32A  -  X.  This  condition  may  be  relaxed 
when  no  density  discontinuity  exists.  The  user  is  encouraged  to  experiment  with  different  mixing 
lengths  and  compare  results  with  other  models. 

b)  The  second  form  of  mixing  function  used  completely  localizes  the  extent  to  within  a 
finite  distance  from  the  interface,  but  is  no  longer  Cx  smooth.  Since  we  must  be  able  to  compute 


the  second  derivative,  it  is  required  to  be  at  least  C2  smooth.  Specifically,  we  choose  a  cubic  spline 


over  the  finite  interval  -Lp  *  C  s  Lp 

In  designing  a  cubic  spline  approximation  for  //(C),  we  have  used  four  sub-intervals  of 
length  Lp/2.  Requiring  continuity  of  the  function  and  its  first  and  second  derivatives,  we  define 


mo 


*»  t  i 

l(1  +  r) 

j  lp 
i  £  i.i.3 

5+z;-3(z;> 

2  C  3 

1 


2  2 


f  ^ 


(2.23) 


The  first  derivative  of  this  function  is  then 


tf'(S)  -  »(W  - 


0 

r<1+r> 

Lp 


1 


l-2(f) 
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f(,-r> 

ap 


-L.s^s  — 


Lp  r  Lp 

"T  ^  T 


(2.24) 


_ 


- 


1 


1 


Note  that  6  ( — — )  -  6(-^-)  -  and  6(0)  -  —  ,  so  it  is  obvious  that 


f*(0<K  -  1  (2.25) 

-00 

and,  therefore, 

lim  6(0  -  6(0  (2.26) 

z,p-*o 
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as  required. 


Finally,  the  second  derivative  is 


4  ,,  L 

4  <• 


H"K>  -  *’(«  -  .  -4<r) 

Ln 


-4(>4) 

4  L” 


P  y.  P 

-  — *5«- 


S*i, 


(2  27) 


Combining  (2.27)  with  (2. 18)  gives  the  formula  for  computing  the  density  potential  function  by 
employing  a  cubic  spline  polynomial  smoothing  function. 


This  function  obviously  has  extrema  at  ^  -  ±-^  given  by  6'  (CCT)  “  ±4  For  the 

Lp 

same  value  of  mixing  length,  the  cubic  spline  polynomial  mixing  function  is  found  to  produce  a 
representation  of  6'  (r)  which  is  narrower  in  depth  and  has  a  larger  amplitude  than  the  hyperbolic 
tangent  mixing  function.  Furthermore,  one  can  easily  show  that  by  defining  this  mixing  length  in 
such  a  way  that  the  spacing  of  the  extrema  for  each  function  agree,  i.e., 

Zfy)  -  (2  x  1.32)Z,ptanh)  -  2.64Z,p,thevaluesoftheextremaare&'^>ly)  »36'^nh). 

Lp 

Thus,  for  the  same  width  between  extrema,  the  hyperbolic  tangent  mixing  function  produces  a 


6'  (z)  representation  with  smaller  extrema  than  the  cubic  spline  polynomial  mixing  function.  This 
is  a  consequence  of  the  infinite  extent  of  the  former  versus  the  finite  extent  of  the  latter  coupled 
with  the  normalization  condition  (2.8). 

As  before,  we  define  a  lower  limit  default  of  Z.  .  -  2A z,  and  when  this  limit  is  reached 

"min 

the  mixing  function  is  centered  on  the  mesh  nearest  the  true  interface  depth.  Applying  the  condi- 
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tion  Lp  &  X  yields  the  condition  Ar,^  ~  - .  This  is  a  stronger  condition  than  before  and  suggests 


that  twice  the  former  transform  size  is  needed  to  utilize  the  polynomial  mixing  function.  However, 
because  of  the  finite  nature  of  this  representation,  this  condition  may  be  relaxed  slightly.  A  good 
rule  of  thumb  for  either  mixing  function  is  Ar,^  -  X.  when  density  discontinuities  are  important 
for  computing  the  correct  reflection  from  the  interface. 

Finally,  let  us  return  for  a  moment  to  the  original  expression  for  Ui(z)  given  by  (2. 16).  If 
we  ignore  the  singular  nature  of  the  second  term  and  assume  that  the  generalized  functions  intro¬ 
duced  are  valid,  we  may  attempt  to  investigate  the  validity  of  the  approximate  form  (2.18)  by 
examining  the  relative  magnitudes  of  the  two  terms  in  (2.16).  When  the  hyperbolic  tangent  func¬ 
tion  is  employed,  the  maximum  value  of  the  first  term,  up  to  a  constant  factor,  has  been  shown  to 


be  approximately  (0.096) 


(P  *-Pj 


Up  to  the  same  factor,  the  second  term  has  a  maximum 


value  of  ^  0ne  can  then  easily  show  that  the  maximum  value  of  the  second  term 

16  (p„+pJV 

is  greater  than  the  maximum  value  of  the  first  term  only  when  pfc  >  3pw.  For  smaller  density  dis¬ 
continuities,  this  analysis  suggests  the  approximate  form  for  Ui(z)  should  be  adequate.  Similarly, 
when  the  cubic  spline  polynomial  is  used,  the  first  term  has  been  shown  to  have  a  maximum  of 

(P<,”PW)  3  (p^  -  pw)  “ 

2 - w  ■  while  the  second  term  has  a  maximum  of  ■=  — - - — r .  A  similar  analysis  shows 

i;  2  (p6+p„)  v 

that  the  second  term  exceeds  the  first  term  only  for  >  5pw  This  implies  a  further  improvement 

in  accuracy  when  the  polynomial  representation  of  the  density  discontinuity  is  selected  instead  of 
the  hyperbolic  tangent. 
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2.2 


Surface  interface 


The  UMPE  model  treats  the  surface  as  a  perfect  reflector  due  to  a  pressure  release  bound¬ 
ary.  This  is  a  Dirichlet  boundary  condition  defined  by 

VU-  0)  -  0  .  (2.28) 

A  popular  technique  used  in  PE  SSF  models  to  achieve  this  is  the  image  ocean  method  developed 
by  Tappert  (1977).  With  this  method,  we  assume  an  identical  image  ocean  overlays  the  real  ocean 
for  negative  values  of  depth  qnd,  furthermore,  the  acoustic  field  is  exactly  equal  but  of  opposite 
sign  in  the  image  ocean,  i.e. 

*(-*)  -  -V(z)  ■  (2.29) 

The  boundary  condition  (2.28)  is  then  satisfied  automatically. 

In  our  numerical  implementation,  therefore,  we  must  define  our  environmental  and  field 
arrays  to  be  twice  as  long  (i.e.  twice  as  deep)  as  necessary  to  describe  the  real  environment  and  real 
acoustic  field.  After  each  range  step,  the  UMPE  model  assures  this  symmetry  by  simply  imposing 
condition  (2.29)  on  the  image  field  fore  <  0.  This  formulation  allows  direct  implementation  of  the 
split-step  Fourier  algorithm  given  by  Eq.  (1 .29)  using  the  full  FFT  transformations  from  e-space  to 
£-space.  It  should  be  noted  that  some  PE/SSF  models  (e  g.,  the  Navy  standard  model)  use  condi¬ 
tion  (2.28)  to  simplify  the  implementation  by  using  only  sine  transforms  (rather  than  cosine 
transforms  since  (2.29)  requires  the  full  field  have  odd  symmetry  about  e  =  0).  This  has  the  effect 

of  reducing  the  necessary  transform  and  array  sizes  by  ^ .  However,  this  simplification  is  only 

valid  when  the  surface  is  completely  flat  and  at  depth  z  =  0.  To  compute  the  exact  forward  scatter 
due  to  a  rough  surface  interface,  as  will  be  discussed  in  a  later  section,  the  full  description  of  the 
field  is  required. 
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3.  Volume  attenuation 


Adding  attenuation  to  the  environment  effectively  introduces  a  complex  component  of  the 
sound  speed  or,  alternatively,  of  the  wavenumber.  The  simplest  way  to  introduce  this  phenomenon 
into  the  wave  equation  is  to  define  the  equivalent  index  of  refraction 


,2  2  . 

n  -  n  +  »a 


(31) 


As  before,  this  can  be  incorporated  into  our  PE  model  by  introducing  a  new  potential  function 


<W*>  - 


.«'(*) 


a(r) 

2  An 


(3.2) 


where  a'  -  —  is  strictly  valid  only  when  the  SPE  approximation  is  invoked.  We  shall  assume  the 

validity  is  more  general,  however,  for  attenuation  values  that  are  not  exceptionally  large.  Based  on 
the  formulation  of  the  algorithm,  this  clearly  has  the  effect  of  damping  the  solution  in  the  space- 
domain  by  the  factor 

e-ikQ&rUlo„(z)  _  ^-Ara'(z)  ^  ^ 

In  terms  of  transmission  loss,  this  reduces  the  field  by 

TLa  -  -20log(e~Ara  )  -  8.686 Ara'dB  .  (3.4) 

The  UMPE  model  assumes  values  input  for  volume  attenuation  have  units  [dB/km/Hz],  Internally, 
these  values  are  rescaled  to  have  units  of  inverse  length,  where  the  unit  of  length  is  specified  by  Ar, 
and  the  field  strength  is  reduced  by  the  factor  given  in  Eq.  (3.3)  every  range  step. 


3.1  Water  volume  attenuation 

There  are  currently  many  common  forms  of  empirical  formulae  defining  volume  attenua¬ 
tion  in  sea  water  as  a  function  of  frequency  and  depth.  A  popular  form  credited  to  Thorp  (1967)  is 
given  by 
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(3.5) 


aT  (dB/kyd)  -  0.003  +  -^4  +  4°^  .  +  2.75  x  10"4/ 

1+/  4100+/ 

where /is  the  frequency  in  kHz.  Transforming  this  to  units  of  [dB/km/Hz]  is  accomplished  by 

a'(r-0)  (3  6) 

The  dependence  with  depth  is  then  given  by 

a'(z)  -  a'  (z  -  0)  [1  -  6.46  x  10-5r(m)  ]  (3.7) 

In  the  UMPE  model,  a  particular  form  for  a'  (r)  is  encoded  and  the  environmental  poten¬ 
tial  function  propagator,  e~,ko^rU°p^  js  scaled  by  the  factor  e“Ara  at  each  range  step.  Currently,  the 
form  of  a'  (z)  is  considered  range-independent.  In  fact,  as  of  the  date  of  this  report,  no  clear  def¬ 
inition  of  a'  (z)  has  been  agreed  upon  and  the  portion  of  the  UMPE  code  computing  a'  (r)  has 
been  commented  out. 

3 2  Bottom  volume  attenuation 

The  UMPE  model  allows  the  user  to  input  a  volume  attenuation  coefficient  for  each  bottom 
profile,  i.e.  a'  may  be  range-dependent  but  is  constant  over  the  depth  of  the  bottom  layer.  Because 
the  volume  attenuation  was  introduced  as  part  of  a  complex  index  of  refraction,  the  discontinuities 
between  the  attenuation  values  over  an  interface  are  smoothed  in  an  identical  manner  to  the  sound 
speed  smoothing  previously  described.  This  is  part  of  the  creation  of  the  array  of  volume  attenua¬ 
tion  values  as  a  function  of  depth,  a'  ( z ) ,  for  all  depths.  As  described  above,  this  is  used  to  scale 
the  values  of  the  environmental  potential  function  propagator. 

3.3  Effective  attenuation  due  to  shear 

As  was  briefly  mentioned  in  a  previous  section,  the  UMPE  model  treats  the  conversion  pro¬ 
cess  of  compressional  to  shear  waves  at  the  bottom  interface  as  a  loss.  The  physical  justification 
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is  that  shear  waves  travel  more  slowly  and  at  steeper  angles  than  the  incident  compressional  waves 
which  generated  them  Consequently,  over  a  given  range  step,  they  attenuate  more  rapidly.  It  is 
unlikely  that  the  shear  waves  will  exist  long  enough  to  be  refracted  or  reflected  back  towards  the 
surface  to  transfer  energy  back  into  compressional  modes.  We  therefore  desire  to  compute  the  frac¬ 
tion  of  incident  energy  which  is  converted  to  shear  wave  energy  at  some  bottom  interface,  and  then 
simply  remove  that  from  the  problem. 

Tapper!  ( 1 985)  derived  an  analytical  expression  valid  for  low  grazing  angle  reflection  at  the 
interface  of  a  lossy  fluid  and  a  lossless  solid  bottom.  The  solid  was  replaced  by  an  equivalent  fluid 
bottom  with  an  effective  attenuation  by  matching  the  reflection  coefficients  of  the  lossless  bottom 
with  and  without  shear.  The  effective  attenuation  was  then  a  function  of  the  bottom  shear  speed 
and  represented  the  simplest  concept  of  loss  due  to  shear  conversion. 

More  recently,  Tindle  and  Zhang  (1992)  developed  a  more  rigorous  approach  for  the  inter¬ 
face  between  a  lossy  fluid  and  a  lossy  solid  bottom  by  attempting  to  match  the  parameters  of  an 
equivalent  lossy  fluid  bottom  by  comparing  the  total  reflection  coefficients  for  each  case.  This  pro¬ 
vides  a  good  approximation  to  both  the  phase  and  amplitude  of  the  reflection  coefficient  for  all 
angles.  They  showed  that,  in  addition  to  an  effective  bottom  attenuation,  it  is  also  necessary  to 
define  an  effective  bottom  density  which  is  smaller  than  the  true  density.  Specifically,  a  bottom 
with  compressional  and  shear  speeds  cb  and  cr  respectively,  density  pb  and  compressional  and 
shear  attenuations  ab  and  a„  respectively,  can  be  represented  by  an  equivalent  fluid  bottom  with 
compressional  speed  effective  density 

P4'  -  Pjfl-2*2^)  ,  (3.8) 

\  (I)  ' 

where  <n  -  2 jcJ  is  the  angular  acoustic  frequency,  and  effective  attenuation 
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cu 


a*  + 


cbcs\K  -  — ; r  •,  1/2  ,  1/2-1 


a»4(l  -2*2^) 
V  co  / 


*1 


; 


(39) 


The  value  to  be  used  for  the  wavenumber  k  remains  somewhat  ambiguous.  For  Eq.  (3  .9), 


U) 


the  authors  chose  k  -  — .  However,  for  Eq.  (3.8)  they  claim  for  several  cases  investigated  it  was 


U) 


found  that  k  -  —  was  the  best  choice  for  when  the  interface  compressional  sound  speed  contrast 
cb 


i/  L 

satisfied  —  <1.2.  For  ratios  exceeding  1 .2,  they  suggest  averaging  the  values  of  effective  density 


CD 


CD 


obtained  by  defining  k  -  —  and  k  -  — .  To  avoid  a  discontinuous  change  in  the  value  of  pb. 


and  because  a  sound  speed  ratio  of  1 .2  is  a  large  contrast  in  most  ocean  environments,  we  shall  only 
employ  the  former  definition.  These  equations  then  simplify  to 


pi'-p‘('-22 


(3.10) 


and 


ab  -  «fc+ 


->  2  1/2  ■)  ,  1/2 
<°(c6-cw)  ( cw~c~s ) 


cbcw 


(3-11) 


Note  that  Eq.  (3.11)  requires  that  cs<cw<cb.  These  formulae  are  implemented  in  the  UMPE 
model  at  each  interface  with  the  understanding  that  the  upper  layer  is  treated  as  the  fluid  and  the 
lower  layer  is  treated  as  the  solid. 
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4.  Rough  interface  forward  scatter 

In  this  section,  we  shall  consider  how  the  UMPE  model  treats  both  rough  bottom  interfaces 
and  surface  displacements.  Because  the  UMPE  model  is  a  one-way  PE  model,  only  forward  scat¬ 
tering  can  be  considered.  Therefore,  it  is  best  suited  for  small  slope  specular  reflection  and  forward 
diffuse  scatter.  Backscatter  from  cliff-tike  structures  can  be  computed  with  the  addition  of  specific 
pieces  of  code  for  individual  problems  but  this  is  not  within  die  general  framework  of  the  current 
version. 

For  both  the  ocean  surface  and  bottom  interfaces,  the  roughness  is  assumed  to  be  charac¬ 
terized  by  a  power  law,  or  fractal,  2-D  spectrum  at  high  wavenumbers  or  small  scales,  i.e. 

H'2(*»  •=-!-)  -  ct*'"  (4.1) 

A corr 

where  Lmrr  is  the  correlation  length  of  the  roughness  and  0  is  the  spectral  exponent.  (Note  that  our 
use  of  wavenumber  is  now  with  respect  to  range  scales  of  interface  roughness.)  To  give  the  full 
spectrum  a  realistically  smooth  structure,  we  assume  the  spectral  form 

HM*) - - - 57;  •  (4.2) 

0*4^) 

The  normalization  factor  p  is  defined  in  terms  of  the  rms  roughness  o  by  requiring 

00 

2  nfW2  (k)  kdk  -  cf2  .  (4.3) 

o 

This  leads  to 

K  -  i(|-Dcr4rr  '  («•<) 

For  the  purpose  of  computing  the  acoustic  field  in  two  dimensions  (depth  and  range),  we 
need  only  the  1-D  roughness  spectrum  along  the  track  of  interest.  Assuming  the  direction  of  prop¬ 
agation  is  along  the  x-axis,  we  need 
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(*,)  -  fw2(Jk;  +  k;)dky 


It  follows  immediately  from  Eq.  (4.3)  that 


^Wx(kx)dkx  -  o2  (4.6) 

-ao 

An  alternative  formula  to  (4.  S)  for  computing  W\  ( k )  along  any  track  may  be  derived  by  converting 


to  cylindrical  coordinates  giving 


WAk)  -  2 1 


W-,  {K)dK 


Substituting  (4.2)  and  (4.4)  into  (4.S)  yields 


-(ill) 

W}(k)  -  Y «2£corrO+l2corS2)  1 


2,P  1W*  dt 

Y’«(2  X)S  ,  p/2 

0(1+0 


(f-Dr(I)r(ti) 


where  T(x)  is  the  gamma  function. 

Within  the  UMPE  model,  we  avoid  computing  gamma  functions  by  assuming  the  1-D  spec¬ 
trum  has  the  form 


»V(*>  -  O+llr/) 


(4.10) 


In  order  to  compute  various  stochastic  realizations  of  roughness,  we  superimpose  upon  this  a  ran¬ 
dom  wavenumber  amplitude  and  phase.  Because  the  complex  amplitude  of  each  wavenumber 

component,  Ae®,  should  exhibit  a  normal  distribution  in  the  complex  A-plane,  the  random  phase 
of  each  component  can  be  obtained  from  ft  -  2nr ,  where  rx  is  a  uniformly  distributed  random 
variable  in  the  interval  (0, 1 ) .  The  magnitude  A,  however,  exhibits  a  Rayleigh  distribution. 
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Because  we  are  creating  a  realization  of  the  power  spectrum,  we  consider  the  magnitude  squared 
which  has  a  negative  exponential  distribution  The  random  amplitude  of  each  component  of  the 

power  spectrum  can  then  be  obtained  by  A"  -  -In  (r2)  where  is  another  independent  uniformly 
distributed  random  variable  in  the  interval  (0,  1 ) .  This  complex  spectrum  is  then  Fourier  trans- 

_  (jc) 

formed  to  yield  a  roughness,  q(.r),  which  is  normalized  by  the  actual  rms  value,  i.e  - — 

<i!  w  > 

The  result  of  this  calculation  is  a  single  stochastic  realization  of  an  interface  roughness  with  the 
proper  spectral  shape  given  by  (4. 10)  and  an  rms  value  normalized  to  unity.  Multiplication  by  the 
desired  rms  value,  o,  should  give  us  exactly  the  rough  interface  defined  by  (4.8)  and  (4.9).  Fur¬ 
thermore,  to  avoid  introducing  interface  displacements  not  found  in  actual  bathymetric  databases, 
the  spectrum  for  bottom  roughness  is  high-pass  filtered  to  remove  wavelengths  larger  than  or  on 
the  order  of  the  database  resolution.  (The  cut-off  for  this  filter  is  specified  within  the  source  code 
and  may  be  edited.)  The  UMPE  model  accepts  as  input  the  total  number  of  realizations  to  compute. 
If  more  than  one  realization  is  requested,  the  calculations  of  transmission  loss  for  given  depths  are 
averaged  incoherently  to  yield  a  measure  of  the  total  rms  field  or  average  power. 

4.1  Water/bottom  interface 

The  treatment  of  the  boundary  condition  at  bottom  or  sediment  interfaces  has  already  been 
discussed  in  detail  in  section  2.  Since  this  treatment  was  independent  of  the  actual  depth  of  the 
interface,  no  additional  calculations  are  necessary  once  the  bathymetry  has  been  specified.  One 
should  note,  however,  that  with  any  type  of  range-dependent  bathymetry,  whether  smooth  sloping 
bottoms  or  with  small  scale  roughness,  each  range  step  of  the  calculation  assumes  a  constant  envi¬ 
ronment.  Recall  that  this  follows  from  the  approximation  J’r(r  +  Ar)  H(r')dr'  •  A rH{r)  implicit 

in  the  PE/SSF  algorithm.  This  does  not  imply  that  the  effects  of  slopes  are  neglected.  If  the  envi¬ 
ronment,  hence  the  bathymetry,  is  sampled  well  enough  and  the  changes  occurring  over  one  range 
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step  are  relatively  small  (e  g.  the  change  in  bathymetry  over  one  range  step,  or  the  slope,  must  be 
less  than  a  wavelength,  6t]  «  XQ)  then  diffusive  effects  will  smear  out  these  apparently  discontin¬ 
uous  jumps  The  combination  of  a  number  of  such  “jumps”  over  several  range  steps  will  be 
observed  to  create  the  anticipated  effects. 

4.2  Surface  interface 

The  calculation  of  rough  surface  forward  scatter  is  more  complicated  because  the  treatment 
of  the  surface  boundary  condition  assumes  that  the  surface  is  fixed  at  z  =  0.  The  UMPE  model  pro¬ 
vides  two  methods  of  computing  rough  surface  scatter.  The  first  maintains  the  imposed  symmetry 
about  z  =  0  by  invoking  an  approximation  which  replaces  the  interface  displacement  by  a  volume 
effect.  The  second  method  is  an  exact  technique  derived  from  first  principles  when  the  surface 
interface  is  not  assumed  to  be  flat. 

4.2.1  Approximate  surface  forward  scatter 

We  assume  the  surface  interface  is  defined  by 

r-rj(r)  -  0  (4.11) 

where  ij(r)  is  the  random  surface  displacement,  increasing  downward,  and  is  defined  as  a  zero- 
mean  Gaussian  random  function  with 

<r,(r))-0  (4.12) 

and 

00 

<tl(r) il (<■')>  -  R ' fWt(k)em'-r)dk  (4.13) 

-00 

This  definition  is  consistent  with  the  interface  roughness  previously  described. 

We  begin  with  the  Helmholtz  equation  for  the  acoustic  pressure  with  a  Dirichlet  boundary 
condition  at  the  surface, 
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and 


V2/?  +  np  -  0 


(4.14) 


/»(n(0)  -  o  .  (4.15) 

It  can  be  shown  that  the  boundary  condition  may  formally  be  put  into  the  volume  as  a  singular 
index  of  refraction.  For  small  surface  slopes,  this  is  approximately 

n2_nW4„(.)1  o  (4.16) 

*0  J 

and  now 

p( 0)  -  0  (4.17) 

Invoking  the  usual  approximations  in  the  standard  parabolic  equation  to  Eq.  (4.16)  leads  to 

+  jrd-A~ko\.u(r'z)  ♦HO'.*)]*  “  0  (418) 

dr  2ko dz~ 

with 

tp(0)-0,  (4.19) 

U{r,z)  is  the  standard  PE  potential  function  (Eq.  (1.13)) 

U(r,z)  - 

A* 

and  p(r,z)  is  given  explicitly  by 

|i(r,.-)  -  l^U'(z)  .  (4.20) 

2k20 

The  variat  le  z  is  now  interpreted  as  the  depth  increasing  downward  from  the  mean  sea  level  atz  = 

0. 

As  before,  we  must  replace  the  functional  6  "(z)  by  some  appropriate  generalized  function. 

The  first  non-vanishing  moment  of  b"  (z)  is  the  second  moment  defined  by 

00 

Jz2b"(z)dz  -  2  .  (4.21) 
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Furthermore,  the  expectation  value  is  given  by 

00  00 

(6"(z))  -  Jy*b"(z)ydz  -  J*6(r)  +  2ip*'tp']</r  -  2  (0)|  (4.22) 

-00  -00 

which  follows  from  an  integration  by  parts.  Introducing  g(z)  as  a  generalized  function,  we  require 
agreement  with  the  expectation  value  of  (4.22),  i.e 

00 

(#(-)}  -  J q*g(z)ydz  -  2  JU|>(0)  (4.23) 

-00 

Since  g(z)  will  have  influence  only  near  the  surface,  we  expand  the  held  in  a  Taylor  series  near  z  = 
0, 

tp(z)  «rJUl>(0)  +...  (4.24) 

To  lowest  order 

00 

(g  (z)  >  -  Jz2g  ( z )  JU|>  ( 0)  |  (4.25) 

-0D 

and  the  condition  on  g(z)  reduces  to  matching  the  second  moments  such  that 


J  z2 g  (z)dz  -  2  . 
-00 

This  can  be  accommodated  by  a  Gaussian  of  the  form 

(  \  A  ~*/L\ 

g(z)  -  Ae 

where  the  amplitude  is  defined  by 

A  fz2e~//L’dz  -  2  . 
-00 

From  a  standard  table  of  integrals  this  reduces  to 


(4-26) 


(4.27) 


(4.28) 


(4.29) 


36 


The  approximate  surface  forward  scatter  is  then  carried  out  by  imposing  the  usual  odd  sym¬ 
metry  at :  =  0  and  adding  to  the  environmental  potential  the  function  p(r,r)  defined  by  (4.20)  with 
the  generalized  function  given  by  (4.27)  and  (4.29).  The  standard  PE  potential  function  was  used 
here  for  simplicity.  However,  the  formulation  of  (4.20)  is  more  general  and  can  be  added  to  what¬ 
ever  environmental  potential  form  is  desired. 

It  should  also  be  mentioned  at  this  point  that  for  surface  interface  roughness  the  UMPE 
model  takes  input  values  of  wind  speed  which  may  be  range-dependent.  An  accepted  empirical 
relationship  (Neumann  and  Pierson,  1966)  is  used  to  determine  the  rms  roughness 


o  - 


(wind  speed)' 

10i 


where  g  ■  acceleration  of  gravity .  The  roughness  correlation  length  is  then  defined  by 


Lcorr  -  200  . 


(4.30) 


(4.31) 


As  before,  the  use  of  a  generalized  function  is  inherently  ambiguous  with  respect  to  die 
choice  of  mixing  length  La  in  (4.27).  Presumably,  this  function  should  not  be  much  wider  than  a 
wavelength  and  yet  be  large  enough  to  allow  adequate  sampling  by  the  given  depth  mesh.  Cur¬ 
rently,  a  minimum  default  value  of  Z,Jmjn  -  A z  has  been  defined.  To  verify  some  range  of  validity 

for  Ls  requires  a  comparison  with  exact  solutions  of  surface  scatter.  A  method  of  computing  such 
solutions  is  the  topic  of  the  following  section. 


4.2.2  Exact  surface  forward  scatter 

We  return  to  the  standard  parabolic  equation  and  impose  the  surface  boundary  condition  at 
the  true  depth  of  the  surface  displacement,  i  .e. 

thl>  1  d2ii> 

‘rr*  -  °  «-32> 

with 

tMr-ri(r))  -  0  .  (4.33) 
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Similar  to  the  previous  treatment  of  extending  the  field  in  an  image  ocean  with  odd  symmetry  about 
z-  0,  we  now  extend  the  field  with  odd  symmetry,  and  U(r,z)  with  even  symmetry,  about  z  =  q(r). 
The  equations  defining  the  field  propagation  in  the  real  and  image  oceans  are  then  (Tappert  and 
Nghiem-Phu,  1985) 

dtp  1  d 

real  ocean:  i-£  +  - j  -*0U(r,z)y>  -  0  ,  (z>q(r))  (4.34a) 

dr  2k0dz 


dil»  chichi)  1  o  il> 

image  ocean:  i(-£  +  2—  ^)  +  f  -k0U(r,z)y  -  0  ,  (z<T)(r))  .  (4.34b) 

It  can  be  shown  that 

(r, -z  +  2t)  (r) )  -  -i|)  (r,z)  =>  i|)  (r,  r\  (r) )  -  0  , 
satisfying  our  required  boundary  condition. 

To  implement  this,  let  us  define  the  total  field  extending  over  the  real  and  image  ocean  as 
ij>  (r,  z) .  Then 


q>(r,z)  ,  (z>q(r)) 

^(r,Z)  "  {  /2A0ti(*-Tl)  /  x  /  ,  xx 

e  °n  1  t|»(r,z)  ,  (z<q(r)) 

.  (r  >  r|  (r) ) 

{_e'2*o^(r-T|)1j,(r>  -z  +  2q(r))  ,  (z<T](r)) 


Using  the  FFT  convention  in  the  UMPE  code, 

ip  (2)  -  FFT  [*(*)]  -  V  . 

and  employing  the  “frequency-shifting”  theorem 

t|)(z  +  z')  -  ^keikiz*z)  -  ^keikxeikz  -  FFT  [fy (k) e**  ] 

Eq.  (4.35)  becomes 


t()(r,z) 


? 


V 


ikz 


,  (z  >  q  (r) ) 


(^~n)  -x  -2ikx\  ikz  ,  /  w 


(4.35) 


(436) 


(4-37) 


(4.38) 
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where  we  have  also  used  the  identity 

*(-*)  -  FFT  [*(-*)]  . 

Eq.  (4.38)  is  implemented  in  the  UMPE  model  to  compute  exact  surface  forward  scatter. 
To  do  so,  an  additional  field  array  has  been  added  to  the  calculation  to  accommodate  the  second 
formula  of  Eq.  (4.38).  At  each  range  step,  before  the  solution  is  transformed  from  A-space  to  z- 
space  a  copy  of  die  field  is  made.  Each  formula  in  (4.38)  is  then  applied  to  the  two  copies  of 

which  produces  two  different  fields  in  r-space.  These  two  fields  are  then  recombined  about  the 
position  of  the  interface  to  produce  the  total  field  prior  to  the  next  range  step. 

Note  that  Eq.  (4.38)  is  also  independent  of  the  type  of  PE  approximation  used  since  both 
(4.34a)  and  (4.34b)  could  have  been  written  in  terms  of  a  general  function  without  affecting 
the  influence  of  die  surface  displacement.  Furthermore,  this  formulation  requires  a  complete  Fou¬ 
rier  transform  as  the  symmetries  involved  cannot  be  accommodated  by  only  a  sine  transform.  A 
drawback  of  this  technique  is  the  necessity  of  additional  FFTs  in  the  SSF  algorithm,  thereby 
increasing  the  total  run  time.  In  this  regard,  the  approximate  forward  surface  scatter  technique  is 
desirable.  However,  to  this  date  no  thorough  investigation  of  the  validity  of  the  approximate  tech¬ 
nique,  presumably  by  comparing  results  from  the  two  methods,  has  been  performed. 

4.3  Effect  of  near-surface  babbles 

The  presence  of  bubbles  in  the  ocean  volume  is  characterized  by  the  probability  density 
function  of  the  number  of  bubbles  per  unit  volume  per  unit  radius,  n  (Jr,  a) .  The  total  number  of 
bubbles  per  unit  volume  is  then 


N{x)  -  Jn(x,a)da  (4.39) 

o 

and  the  void  fraction,  defined  as  the  volume  of  bubbles  per  unit  volume,  is 
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(4.40) 


V(x)  -  ^-xan(x, a)da  . 
o 

We  shall  assume  the  primary  effect  of  the  presence  of  bubbles  is  to  alter  the  acoustic  index 
of  refraction  which  is  now  defined  by 

«(*)  -  »w($)  +»*(*)  (4  41) 
nw  (x)  is  the  usual  index  of  refraction  of  water,  and  nb  (x)  is  the  effective  index  of  refraction  due 
to  bubbles  within  the  water  volume  defined  by 

nb(i)  -  const  x  V(x)  (442) 

where  the  proportionality  factor  is  some,  as  yet,  undetermined  constant.  Presumably,  the  bubbles 
are  organized  in  such  a  way  that  the  average  void  fraction  is  a  function  of  depth  only,  i.e. 

V(x)  -  V(z)  +6F(x)  (4.43) 

and  we  shall  assume  that  the  bubbles,  hence  nb  (x) ,  are  concentrated  near  the  surface.  We  may 
then  expand 

nb(*)  mnbo(x>y)&(2)  +nbi(x,y)b'(z)  +nb2(x,y)bH (z)  +  ...  (4.44) 

When  the  PE/SSF  algorithm  is  applied  to  this  expansion,  the  odd  terms  vanish  since  n(x) 
must  be  an  even  function  about  z  ~  0  (and  the  odd  derivatives  of  the  delta  function  are  odd  func¬ 
tions).  The  first  two  even  terms  in  (4.44)  are  defined  by 

00  00 

»bo(x,y)  -  Jnb(x,y,z)dz  -  2jnb(x,y,z)dz  (4.45) 

-oo  0 

and 

00  00  00 

nb2  (*,y)  -  j  J znb  (x,  y,  z)  dz  -  fz2nb(x,ytz)dz  -  const  x  Jz2V(x,y,  z)dz  (4.46) 
-00  o  0 

The  first  term,  however,  will  also  vanish  since  the  delta  function  combined  with  the  Dirichlet 

boundary  condition  at  z  *  0  has  no  effect.  The  leading  term  in  the  correction  for  the  presence  of 

bubbles  is  then  nb2{x,y)bn{z). 
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We  now  recall  in  our  development  of  the  approximate  surface  forward  scatter  there  was 
also  a  term  proportional  to  6"  (z)  defined  by  Eq.  (4.20).  It  therefore  seems  plausible  that  we  may 
incorporate  the  effects  of  near  surface  bubbles  into  our  forward  scattering  calculations  by  defining 
an  effective  surface  displacement, 

00 

tleff  (*,y)  -  r\  (x, y)  +  const  x  AqJ z2  V (x, y,  z)  dz  .  (4.47) 

o 

Using  this  approach,  we  are  then  free  to  employ  either  forward  surface  scatter  method 

There  is  presently  no  clear  definition  for  the  constant  proportionality  factor  in  Eq.  (4.47). 
Furthermore,  there  is  no  accepted  empirical  relationship  between  wind  speed  and  bubble  void  frac¬ 
tion.  It  is  not  clear  how  we  should  define  these  quantities  in  terms  of  the  UMPE  model  input 
parameters.  Currently,  the  model  allows  the  user  to  input  simply  an  effective  rms  roughness  and 
correlation  length  characterizing  the  virtual  surface  displacement  due  to  bubbles.  With  these  val¬ 
ues,  a  “bubble  displacement”  is  computed  using  a  spectrum  of  the  same  form  as  the  rough  surface 
but  with  independent  amplitude  and  phase  realizations.  The  total  surface  roughness  is  then  formed 
by  combining  the  surface  and  bubble  displacements  together. 

Incidentally,  it  is  also  unclear  whether  each  random  interface  realization  should  be  indepen¬ 
dent.  Specifically,  should  surface  waves  and  near-surface  bubbles  have  similar  spectral  shapes  and/ 
or  phases?  One  may  also  consider  the  correlation  between  interfaces  of  a  thin  sediment  layer  and 
the  underlying  basement.  These  possibilities  are  currently  ignored. 


5.  Broadband  travel  time  analysis 

The  calculation  of  the  time  domain  arrival  structure  at  a  given  range  can  be  performed  in  a 
straightforward  manner  using  a  PE  code.  Recall  that  we  began  this  report  by  assuming  a  form  for 
the  time-harmonic  acoustic  field  defined  by 

P  (r,  z,  co/)  -  pm  (r,  z)  e~,mt  .  (5. 1) 
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The  representation  of  the  field  in  the  time  domain  is  then  simply 

P(r,z,t)  -  FFT[>w(r,r)]  -  ^/>w(r, z)  e~'m  (5.2) 

which  follows  by  recognizing  that  (5. 1)  is  a  single  (CW)  component  of  the  general  time  dependent 
field  (5.2).  In  other  words,  by  choosing  a  particular  frequency,  id',  we  are  computing 

-  ^Pm(r,z)e~m,bOKO),  -  P(r,zt  ta’t)  (5.3) 

at  time  /  *  0. 

To  compute  the  arrival  time  structure  at  some  fixed  range  r  =  A,  we  need  to  compute  the 


complex  field  p^xiz)  for  many  frequencies  and  then  Fourier  transform  to  obtain  PK  t) ,  die 
set  of  complex  pressure  values  in  time/depth  space.  Because  the  FFT  assumes  inputs  over  the  fre- 
BW  BW 

quency  band  f0  — to  /0  +  — ,  where  EW  is  the  bandwidth  of  the  acoustic  source  and  /0  is 

die  center  frequency,  it  becomes  computationally  burdensome  to  consider  high  frequency  calcula¬ 
tions  even  if  the  frequency  bandwidth  of  computed  fields  is  small.  Therefore,  the  UMPE  model 
heterodynes  the  signal,  shifting  the  center  frequency  to  d.c  In  other  words,  we  compute 

>*<-•.')  -  (5<) 

where  (d0  -  2j ifQ.  This  allows  us  to  place  the  complex  pressure  values  into  frequency 

bins  symmetrically  about  o>  -  <dq  -  0.  The  inputs  then  extend  over  the  band- ^BW  to+^BW. 

This  has  no  effect  on  the  computed  intensities  which  is  evident  by  noting  that  /  -  P*P  is  inde¬ 
pendent  of  (Do- 


In  practice,  we  compute  the  FFT  of  -p 

va 


w  R  (z) .  Because  this  neglects  the  overall  phase 


.  R 

iknR  'V 

factor  e  -  e  °,  die  time  domain  is  heterodyned  around  the  value  t0 


— .  Arrival  times  are 

co 


then  given  as  values  of  “reduced  time”  or  (/  - 10)  For  time  domain  analysis,  the  UMPE  model 


requires  the  user  to  input  a  center  frequency,  a  bandwidth,  and  the  number  of  frequency  bins  to 
compute  (the  frequency  transform  size).  This  obviously  introduces  some  ambiguity  into  the  abso¬ 
lute  travel  time  prediction  via  the  choice  of  reference  sound  speed.  It  is  expected  that  an 
implementation  of  the  Co-insensitive  model  would  remove  this  ambiguity. 

To  this  point,  we  have  assumed  each  frequency  component  of  the  total  field  contributes 
equally  in  amplitude.  This  is  synonymous  with  assuming  the  frequency  spectrum  of  the  source  is 
flat  over  some  frequency  band,  i.e.  a  “box-car”  function.  It  is  well  known  that  the  transform  of  such 
a  function  will  introduce  side-lobes.  Furthermore,  since  the  transform  size,  hence  the  time-win¬ 
dow,  must  be  finite,  wrap-around  or  aliasing  may  occur.  To  minimize  these  effects,  the  model 
creates  a  sine-squared  filter  which  is  applied  to  the  upper  and  lower  1/8  of  the  bandwidth.  This 
method  of  tapering  greatly  reduces  the  effect  of  side-lobes.  (More  about  similar  filters  used  in  the 
model  will  be  discussed  in  a  later  section.) 

The  user  should  still  be  wary  of  wrap-around  effects.  These  may  arise  if  the  time  window 
is  not  wide  enough.  To  avoid  this,  the  user  should  estimate  the  anticipated  time  spread  of  the  signal 
before  generating  a  solution  and  input  the  proper  parameters  to  allow  for  this  amount.  The  total 
length  of  the  time  window  is  given  by 


T  - 


N 

BW 


(5.5) 


where  BW  is  the  bandwidth  of  the  source  to  be  modeled  and  N  is  the  number  of  frequencies  for 
which  the  model  should  compute  solutions  to  the  field,  or  the  frequency  FFT  size.  The  frequency 
resolution  is  obviously  given  by 


A  co  - 


BW  1_ 

N  "  T  ‘ 


(5.6) 


Thus,  the  wider  the  time  window  needed,  the  finer  the  bandwidth  must  be  resolved  and  the  more 
total  runs  at  different  frequencies  must  be  performed.  This  can  obviously  lead  to  extremely  long 
computer  run-time.  Careful  considerations  should  therefore  be  given  before  entering  these  param 
eters. 
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It  should  be  noted  that  while  the  model  creates  a  rather  simple  tapered  box-car  function  to 
represent  the  frequency  spectrum  of  the  source,  the  code  could  easily  be  modified  to  accept  a  spe¬ 
cific  frequency  spectrum  defined  by  the  user.  Because  the  output  is  in  dB  units  of  transmission  loss, 
all  that  would  be  required  is  a  source  spectrum  normalized  to  maximum  values  of  unity  and  sam¬ 
pled  at  the  appropriate  frequencies.  This  function  would  then  serve  as  the  filter  function  in  the 
frequency  domain.  Furthermore,  the  model  is  unable  to  treat  multiple  interface  roughness  broad¬ 
band  calculations.  If  one  desires  to  compute  an  ensemble  average  of  the  travel  times,  numerous 
broadband  calculations  for  single  stochastic  realizations  could  be  performed  and  the  results  com¬ 
bined  in  a  post-processing  routine. 

Finally,  there  is  a  built-in  ambiguity  in  the  computed  levels  of  transmission  loss,  or  TL  (to 
be  defined  later),  in  the  time  domain.  This  is  due  in  part  to  variability  with  respect  to  the  frequency/ 
time  transform  size,  beamwidth  of  die  source  function,  type  of  source  function  used,  and  several 
other  parameters.  To  ob^in  true  TL  levels,  we  can  simply  compute  the  travel  time  values  out  to 
the  reference  range  of  R0  -  1  m .  (More  about  this  parameter  in  a  later  section.)  The  peak  value 
at  this  range  should  be  TL  -  0  dB  re  lm.  Subtracting  the  actual  value  obtained  from  subsequent 
calculations  should  provide  true  TL  levels.  Note  that  this  is  only  true  when  all  source  and  broad¬ 
band  parameters  remain  the  same. 

6.  Acoustic  particle  velocity 

The  ability  to  compute  acoustic  particle  velocities  from  the  pressure  field  has  recently  been 
added  as  an  option  in  the  UMPE  model.  The  motivation  was  simply  a  lack  of  available  fully  range- 
dependent  acoustic  propagation  models  which  can  produce  this  result.  The  implementation  is 
straightforward  but,  as  will  become  clear,  adds  a  non-negligible  amount  of  run-time  to  perform  the 
necessary  calculations. 

We  begin  by  stating  the  general  form  of  the  parabolic  equation. 
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th f 

Tr  *  ’'VV»  • 


(6.1) 


where  if  is  related  to  the  acoustic  pressure  by 


»  p  r*<W*°r 


(6.2) 


The  operator  can  take  any  of  the  five  forms  defined  in  section  1 .  By  conservation  of  linear 
momentum,  the  acoustic  particle  velocity  $  -  (v^  vz)  is  related  to  the  pressure  by 


v - '—Vp, 

<“Po 


(63) 


assuming  both  v  and  p  have  the  harmonic  time  dependence  e  .  Separating  the  components  of 

v  leads  to  the  expressions 

/  dp  _ 

Vr  "  _U)P0dr 


«*Po 


po  P^o  j? 
CqPoW*  ,2V 


+  4 


and 


(6.4) 


_ _ ‘-*E  „  -  .  (6.5) 

2  Oip0dz  c0p0*0^p0r  dz 

In  the  r-direction,  it  is  found  that  the  application  of  die  operator  must  be  applied  to  the 

field  function  if  in  order  to  perform  the  calculation.  As  previously  noted,  each  of  our  PE  options 
separates  the  operator  Hop  into  two  parts, 

Hop  -  V-^>  *  ■  (6  6) 

dr4- 

Additionally,  the  vertical  component  of  acoustic  particle  velocity  requires  application  of  the  oper¬ 
ator  ^  on  the  if  field.  To  perform  these  operations,  we  consider  the  wavenumber  representation 
of  the  field  at  some  arbitrary  depth  z  defined  by 
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tp(-)  -  FFT  [*(*)]  - 


(6  7) 


then 


fit  (-)  -  -  FFT  [i*i  (*)  ]  (6.8) 

and 


(-)  -  -y^V'**  -  FFT  [-t2 (i)  ] 


(6.9) 


We  now  state  explicitly  the  various  formulas  for  acoustic  particle  velocity  corresponding  to 
the  type  of  PE  approximation  chosen: 


*)  SPE:  r„,  -  -Ag-ij ,  t/.,.-i(»‘-l). 


1  /  J2 

2 


so 


vr(r,z)  - 


coPo 


+  +  ;  (610) 


1*. 


2)  VSPE:  Top 


±j£ 

2*o  dz2 


-(«-  1)  , 


so 


,  x  PQ  K  <*(/ 
vr(r,z)  -  -  - e 

coPo\  Por 


'Tfr+'n‘’(r'-')  -FTT  (5^  (<•.*)] 

l2  *5  ) 


(611) 


2f  2  1/2  l"1 

3)  VWAPE:  Top.  ^(1+^)  +  1 1  -i  (rt2  -  I)  . 


SO 
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vr{r,z) 


P0  f>*0  ik0r\  .y(r,z)  \  2 

-  e  i— =-r +  x(« 

C0P0 1  Por  2V  “ 


^-  +  i(W2+l)tt.(r,z)-FFT(  1-  1-^  $(r,k))  ;  (6 

o'  -  V  \  ki  ) 


4)  WAPE:  T  -  (1  +  -^)  ♦  1  ,  U  -  -(»-  D 

dr  L  dz  J 


•Vfo*)  “  ro~J^e'k°r  '^Tkls  +ttV(r’2')  -PF7  (  1  -  I1  -  ^  ;  (6.13) 

coPoNPor  [  2*°^  V|_  V  *oJ 


5)  LOGPE:  r„,  -  -/Jcm  (-£■£)]  ,  Uop  -  -/»(»)  , 

SO 

vr(r,z)  - 

rrr  +  (!  +/«(w))t|>(r,z)  +FFT(/»rcospli|>(r,*))l  (6.14) 

coPo\Por  L  2*or  L  *oJ  J 

For  any  of  the  above  options,  the  r-component  of  acoustic  particle  velocity  is  computed  from 

v2(^-)  -  •  (6i5> 

coP0SPor  L*o  J 

It  is  obvious  from  these  expressions  that  for  every  range  step,  additional  FFT’s  must  be  computed 
thereby  increasing  the  overall  run-time  by  a  significant  amount. 
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7.  Source  functions 


In  this  section  we  shall  define  the  initial  conditions  for  the  PE  field  function,  tp(r  m  0,  r) 
Previously,  we  have  assumed  the  relationship  between  and  the  acoustic  pressure,  ignoring  the 
effect  of  density,  is  of  the  form 


We  obviously  cannot  define p  at  range  r-  0  from  this  exptession.  This  is  analogous  to  the  unde¬ 
fined  amplitude  of  a  point  source  Green's  function  at  die  source  location.  Therefore,  we  choose  to 
define  the  source  amplitude  relative  to  that  at  some  small  but  finite  distance  from  the  source.  Spe¬ 
cifically,  we  choose 

p(r-R0)  -  P0  (7.2) 

Consistent  with  reference  values  used  in  most  sonar  equations,  we  define  the  reference  range 

/?0  -  1  m  (7.3) 

and  the  source  level,  SL,  is  related  to  PQ  by 

SL  -  20log^  dBre^0  (7.4) 

The  dB  units  of  SL  are  explicitly  stated  relative  to  a  reference  pressure  value  of  Pr  -  1  pPa  at  the 
reference  range  Rq. 

We  are  still  left  with  the  task  of  determining  a  form  for  the  source  field  (r  -  0,  z) .  We 

begin  by  writing  (7. 1)  as 

ii>(r,r)«?'*°r  -  2-JjLp(r,z)  (7.5) 

from  which  it  follows  that 

tj)(r  -  0,  z)  -  lim i-^-  ffp(r,z)  (7.6) 

r-*0P0^KQ 
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In  the  vicinity  of  a  point  source,  we  know  the  pressure  field  takes  the  form  of  the  spherical  Green’s 
function.  Since  tp(r^r)  is  only  in  two  dimensions,  we  write 

P  •  gge*** ,  R  -  JS77 ,  a  -  P0R„  (TT) 

/>„ 

where  a  is  defined  by  requiring  \p\  -  —  at  R  -  R0.  We  represent  the  source  at  ( 0,  zs)  as  a  point 


source  by  defining 


g)(r  -  0,2)  -  a6  (z-rs)  . 


Integrating  both  sides  of  (7.6)  over  all  depths  yields 


a  -  a  lim  f  -^—zek</dz 

r-*0PQ  y Rq  J  2 JlR 


Because  we  are  ultimately  interested  in  the  solution  in  the  far-field  (r  »  2),  we  may  approximate 


R  -  Jr2 +  z2~r+?r-,  2  -  2 - 2S  . 


(7.10) 


Eq.  (7.9)  then  reduces  to 


a-^lim  I--  fe  °2'dz  -  J^lim  i-,  ~ 
*  °r- 0  2  n  Jr  J  V  °r-*0  2xJ?‘l  '*0 


1 2x*o  ' 


(7.11) 


It  is  desirable  to  begin  the  calculation  by  specifying  the  source  in  the  ^-domain  Including 
the  influence  of  tl.  mage  source,  a  straightforward  Fourier  transform  of  (7.8)  yields 


i|»(r-0,  k)  -  a  J*  [&(--£$)  -  6  (2  +  zs)  ]  e~ikzdz  -  -2/asin  (kzs)  .  (7. 12) 


which  indicates  that  the  wavenumber  representation  of  the  starting  field  has  a  constant  amplitude 
modulated  by  a  phase  due  to  the  interaction  of  the  source  and  its  image.  This  constant  amplitude 


is  consistent  with  the  notion  of  an  omnidirectional  point  source  which  puts  equal  energy  into  all 
wavenumbers  (i.e.,  all  directions). 

As  an  approximation  to  the  delta  function  starting  field,  we  shall  consider  a  Gaussian  func¬ 
tion  defined  by 


2  /  2 


/A  x  Cl  -(Z  -I*)'/* 

V  (0,r)  -  — re 

wjn 


(7.13) 


where  the  additional  factor  is  required  to  give  the  correct  normalization  in  the  limit  w  -*  0 .  Includ¬ 
ing  the  image  field  and  transforming  this  to  the  ^-domain  yields 


ty(0,  k)  -  — = 


a 

wja 


S' 


w2  -ik(z  +  zs) 


dz- 


-2 

2 

i  w  “ 

e  e 


ik  (z  -  zs)  j- 


*V 


-  -2iae  4  sin(Jtrs)  .  (7.14) 


The  first  thing  we  notice  about  this  expression  is  that  the  Gaussian  normalization  factor  has  been 
removed  and  the  remaining  factor  is  identical  to  that  of  Eq.  (7. 12)  including  the  phase  modulation 
due  to  the  source/image  interference.  In  this  approximation,  however,  the  constant  amplitude  of 
(7. 12)  has  been  replaced  by  a  Gaussian  of  width  (2/w)  and  maximum  value  at  k  -  0  of  unity.  This 
may  be  interpreted  as  placing  the  correct  amount  of  energy  at  k  =  0  (corresponding  to  horizontal  or 
9  =  0  propagation)  and  quickly  tapering  off  fire  energy  placed  at  higher  absolute  wavenumbers. 
This  is  consistent  with  the  small  angle  validity  of  the  standard  parabolic  equation.  Unfortunately, 
for  angles  significantly  removed  from  horizontal,  this  type  of  source  function  creates  a  poor  repre¬ 
sentation  of  an  omnidirectional  source. 

For  the  wide  angle  PE  approximations,  it  would  seem  reasonable  to  simply  replace  the 
Gaussian  function  by  unity  for  all  wavenumbers  thereby  equally  populating  all  directions  of  prop¬ 
agation.  However,  even  the  wide  angle  approximations  are  assumed  valid  only  up  to  angles  of  40° 
or  so.  Additionally,  the  finite  FFT  size  will  restrict  how  laige  k/kQ  can  be.  Therefore,  a  smooth 
taper  is  included  at  high  absolute  wavenumber  values  to  limit  the  angular  width  of  the  source  func¬ 
tion  and  to  reduce  the  influence  of  sidelobes.  Thomson  and  Bohun  (1988)  have  also  shown  that  a 
wide  angle  source  needs  to  be  modified  by  the  factor 
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(7.15) 


-1/4 


1  ■”*  '  n 

F(k)  -  ( 1-0  .  (1*1  <*0) 


to  produce  the  correct  solution  in  the  far-field.  Note  that  k  -  k0  corresponds  to  0  -  90°,  so 
\k\  >  k0  represents  imaginary  angles  of  propagation  (evanescent  modes).  It  is  required  then  that 


the  source  function  be  tapered  within  the  limits  of 


<  1. 


Finally,  we  consider  the  starting  field  appropriate  for  a  line  array  of  sources.  Assuming  the 
array  can  be  approximated  by  a  continuous  line  array  of  length  L  (the  separation  between  adjacent 
sources  must  be  «  XQ),  it  is  easy  to  show  the  directivity  of  such  a  source  is  defined  by 


D  - 


sin  ( ^-w) 
*o 
k 

-r-W 


W  - 


k0L 


(7.16) 


As  required  from  our  previous  comments,  this  has  a  maximum  value  of  unity  at  * = 0  (horizontal). 
However,  the  amplitude  of  this  type  of  source  function  must  be  altered  to  account  for  the  direc¬ 
tional  enhancement.  In  other  words,  the  source  amplitude  should  correspond  to  an  equivalent  point 
source,  so  we  must  modify  the  source  level  by  the  directivity  factor  corresponding  to  a  linear  array, 

D1  -  101og£>  (7.17) 

where 


by 


dQ 


(7.18) 


The  source  level  of  an  equivalent  point  source  is  then  SL'  -  SL  +  DI,  so  the  amplitude  factor  of 
the  source  function  becomes 

,1/2 


a'  -  aQl 


(7.19) 


Assuming  k0L  »  1 ,  Eq.  (7.18)  yields 


^  2 L  *oL 

q~t-  t  • 


(7.20) 
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The  three  types  of  sources  described  above  are  available  options  for  generating  a  starting 
field  in  the  UMPE  model.  Before  continuing,  we  shall  review  these  functions  which  shall  now  be 
written  in  the  general  form 

tM*)  -  “/(£>  (7.21) 

*0 

The  subscript  “0”  on  the  A-domain  starting  field  is  meant  to  indicate  that  this  formulation  applies 
only  for  a  free-field  source  at  depth  zs  -  0.  We  shall  consider  the  effects  of  zs  0  0  and  the  image 
source  shortly.  Furthermore,  we  now  demand  that  both  a  and  f(k  /  Aq)  are  pure  real  quantities  (such 


that  a  ■  |a|).  This  implies  even  symmetry,  (-k)  -  yQ(k) .  The  three  starting  field  options 
can  now  be  written  explicitly  as 

1)  Gaussian  source: 


/(x)  -  e'*'"' 


a  -  A ■ 
2\*k0)  ’ 


(w  m  user  input  width ) , 

1/2 


(7.22) 


2)  wide  angle  source: 


/(*)  -  (1-x") 

n  1  /2 

-  *4 

2\xk0) 


-1/4 


(7.23) 


3)  line  array  source: 


fix)  - 


sin  (xw) 


xw 
1/2 


w  - 


k0L 


1  / 2Rq\  “  k0L 


1/2 


(7.24) 
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Note  that  an  additional  factor  of  AA  has  been  included  in  the  normalization  This  is  needed  for  the 
discrete  Fourier  transform  which  assumes  both  functions  ty(k)  and  ip  (z)  are  unitless. 

We  now  consider  the  general  form  for  the  complete  Jt-domain  starting  field.  If  we  allow 
“steering”  of  the  starting  field,  such  that  the  /t-domain  formulation  is  about  some  wavenumber 
k  -  kc  -  A:0sin8(,»‘ 0,  then  it  is  fairly  straightforward  to  show  that  the  A-domain  starting  field  for 
a  source  at  depth  z  =  zs  is  given  by 


*<*)  -e"N>0(*-*c) -«%><*+*,)  •  (725) 

It  is  apparent  from  this  formula  that  the  complete  ^-domain  field  has  odd  symmetry, 
ip  (-it)  -  -ip  (k) ,  and  therefore  the  r-domain  field  must  also  have  odd  symmetry  about  z  -  0  as 
required  by  the  surface  boundary  condition.  Eq.  (7.25)  together  with  Eqs.  (7.21)  through  (7.24) 
constitute  the  complete  formulation  of  die  available  starting  fields  in  the  UMPE  model. 

One  may  note  an  apparent  ambiguity  in  die  prescribed  implementation  of  a  steered  source 

when  the  wide-angle  source  is  used  since  ip  (k-kc)  -  ip  (k  +  kc)  -  1  (neglecting  the  wide- 
angle  correction  factor  (7. 1 5)).  This  ambiguity  arises  when  we  attempt  to  steer  a  source  which  rep¬ 
resents  an  omnidirectional  source.  We  have  removed  this  ambiguity  by  allowing  the  source  filter 
function  to  be  steered,  thereby  centering  the  filter  about  k  -  kc. 

As  a  final  note,  we  must  address  the  issue  of  mesh  point  symmetry  in  the  UMPE  model. 
The  model  assumes  the  usual  mesh  symmetry  in  the  ^-domain. 


*(/) 


(/'  —  1 )  AA ,  1  s  i  s  ^  +  1 
1)A*,  ^  +  2  sisAf  . 

A* 


(7.26) 


However,  in  the  z-domain  the  UMPE  model  employs  what  is  sometimes  referred  to  as  the  “half¬ 
mesh”  symmetry  defined  by 


~(0 


Az,  lsis- 

1- (A-  /  +  ^)  A:  ,  J+ls/siV  . 


(7.27) 
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Another  way  to  contrast  these  symmetries  is  by  noting  that  ty(k(i  -  2) )  -  -$>  ( k  ( i  -  AO ) , 
ty(k(i  -  3) )  -  -ip  (£(/-  A  -  1) )  .etc.  while  in  the  z-domain,  tp  (z (/  =  1))  =  -ip  (z(/  =  AO), 
ip  (z(/  -  2) )  -  -ip  (z  ( /  -  N -  1 ) ) ,  etc.  To  determine  how  this  affects  the  source  function,  con¬ 
sider  the  following.  If  we  computed  the  source  function  according  to  (7.25),  the  transform  of  this 
function  would  produce  a r-space  starting  field  ip  (z)  with  the  same  “whole-mesh”  symmetry. 

Instead  of  ip  ( z ) ,  we  wish  to  compute  ip  (z  +  i  Az) .  Therefore,  we  need  to  make  the  replacement 

ik— 

*<*>-*(1/2)  <*>  -*<*)«  2  •  (7.28) 

Our  final  form  for  the  Jt-domain  starting  field  is  then  Eq.  (7.2S)  corrected  for  the  r -domain  symme¬ 
try  by  (7.28).  Note  that  this  is  the  only  time  such  a  “correction”  needs  to  be  made.  Once  this 
symmetry  is  established,  it  remains  throughout  the  calculation.  Of  course,  to  be  consistent  the  z- 
domain  potential  function,  hence  the  environment  and  all  other  calculations  in  the  z-domain,  must 
be  computed  using  the  half-mesh  symmetry. 

In  addition  to  specifying  one  of  the  analytical  formulas  for  the  starting  field  listed  above, 
the  user  may  also  input  a  generic  source  function.  The  input  format  should  be  compared  with  the 
source  code’s  expected  format.  Also,  the  model  will  assume  the  correct  half-mesh  symmetry 
already  exists,  and  the  source  has  been  normalized  to  yield  a  unit  source  level  at  R 0.  (The  true 
source  level  is  already  an  input  parameter.) 

8.  Filters  and  sponges 

As  has  been  mentioned  briefly  in  previous  sections,  there  are  several  places  within  the 
UMPE  model  which  require  the  use  of  a  filter,  or  “sponge,”  to  produce  the  necessary  removal  of 
energy.  The  function  used  in  the  model  is  a  sine-squared  function  which  is  smooth,  has  a  contin¬ 
uous  derivative,  and  goes  from  zero  to  unity  within  a  finite  region.  The  most  obvious  need  for  a 
filter  we  first  recognize  is  the  radiation  condition  ip  (z)  -*  0  as  z  -*  ±00 .  Because  the  computa- 
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tional  depth  is  finite,  however,  we  must  force  the  field  amplitude  to  zero  at  the  maximum  depth. 
Note  this  also  serves  to  eliminate  wrap-around  between  the  real  and  image  oceans.  This  must  be 
done  in  a  smooth  and  relatively  “slow”  manner  to  avoid  reflections  from  the  filter  function  itself. 
As  is  common  practice  in  many  PE/SSF  models,  the  UMPE  model  applies  the  sine-squared  filter 
function  to  the  bottom  quarter  of  the  computational  depth  of  the  real  ocean.  The  user  is  responsible 
for  ensuring  that  this  does  not  interfere  with  propagating  energy  that  would  affect  down-range 
results,  i.e.  the  bottom  depth  should  be  well  above  the  region  of  the  filter.  The  net  effect  of  this 
bottom  filter,  or  bottom  sponge,  is  to  remove  energy  which  has  penetrated  deep  within  the  bottom 
layer  before  it  reaches  the  maximum  computational  depth. 

The  bottom  sponge  is  a  necessary  part  of  the  calculation  to  match  the  radiation  condition 
and  is,  therefore,  always  computed.  Occasionally,  one  may  wish  to  focus  attention  on  bottom  and/ 
or  sediment  interface  interactions.  To  isolate  these  phenomena  from  interfering  surface  reflected 
energy,  it  is  desirable  to  “remove”  the  surface.  The  UMPE  model  allows  this  option  and  accom¬ 
plishes  the  effect  by  copying  the  bottom  sponge,  inverting  it,  and  placing  it  at  the  surface  interface. 
The  result  is  a  sine-squared  filter  in  the  upper  quarter  of  the  real  ocean  depth.  Again,  the  user  is 
responsible  for  ensuring  that  the  region  of  interest  lies  between  the  surface  and  bottom  sponge 
regions.  Once  these  sponges  have  been  defined,  they  are  included  in  the  calculation  of  the  z-space 

propagator  function,  e“'*°ArtV(z)  35  a  multiplying  factor. 

In  addition  to  requiring  t|>  (r)  -*  0  as  z  -*  ,  we  must  also  require  %  (A)  -*  0  as 

k  -*  ±Amax .  The  exact  same  filter  used  to  create  the  bottom  sponge  is  applied  to  the  ^-domain 

propagator  function,  e~,k°*rT°p(~k'>  This  results  in  a  sine-squared  filtering  of  the  upper  quarter  of 
the  positive  k  wavenumber  spectrum .  This  filter  is  then  inverted  symmetrically  about  k  =  0  and 
applied  again  to  produce  filtering  of  the  lower  quarter  of  the  negative  k  wavenumber  spectrum.  In 
this  way,  wrap-around  in  fc-space  is  avoided. 

Finally,  we  shall  consider  the  filter  applied  in  the  creation  of  the  source  function.  In  addi¬ 
tion  to  the  center  angle  and  angular  width,  the  input  includes  the  angular  tapering  width  of  the 
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source  function.  With  these  parameters,  a  filter  function  is  computed  which  has  a  value  of  unity 
over  the  wavenumber  range  corresponding  to  the  given  angular  width  and  centered  at  the  wave- 
number  corresponding  to  the  center  angle  of  the  source .  At  each  end  of  this  unity  function  is  placed 
a  sine-squared  function  of  wavenumber  width  corresponding  to  the  given  angular  tapering  width. 
This  filter  is  applied  to  the  final  form  of  the  source  function  as  described  in  the  previous  section. 
Note  that  this  filter  has  little  effect  if  a  Gaussian  source  function  is  specified  since  the  source  ampli¬ 
tude  will  already  be  reduced  by  e~X  before  the  taper  begins  to  apply.  Similarly,  an  array  source  is 
a  narrow  wavenumber  source  and  is  not  significantly  effected  by  this  filter.  The  main  role  of  the 
source  filter  is  to  force  %  (A)  -*  0  as  A  -*  for  the  wide  angle  source.  It  also  allows  the  wide- 

angle  source  to  be  steered,  creating  greater  versatility  for  this  source  option. 


9.  Calculations  of  transmission  loss 


The  calculation  of  transmission  loss  has  previously  been  defined  by 


TL  (r,  z)  -  lOlog  -  lOlog 

\n) 


PR0  2 

L^1^1 


(9.1) 


Since  the  PE/SSF  algorithm  generates  tJ>(z)  sampled  every  A z,  the  values  for  TL  are  also  sampled 
every  A s  at  each  range  step  separated  by  A r.  This  grid  of  TL  values  can  be  output  to  allow  the  user 
to  view  the  entire  TL  field  which  corresponds  to  the  total  pressure  field.  This  requires  an  external 
plotting  program  be  used. 

In  those  situations  where  we  desire  to  determine  the  transmission  loss  at  a  fixed  depth  (for 
example,  at  a  fixed  receiver  depth  or  at  the  sediment  or  bottom  interface)  we  must  somehow  inter¬ 
polate  the  gridded  TL  values  to  this  depth.  This  is  performed  using  a  five-point  polynomial 
interpolation  of  the  ty(z)  field.  The  UMPE  model  allows  the  user  to  request  the  TL  values  output 
at  either  bottom  interface  and  at  up  to  ten  arbitrary  depths.  Since  a  five-point  interpolation  is  to  be 
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used,  it  is  required  that  all  depths  be  deeper  than  Actually,  depths  as  shallow  as  ^Az  are 


allowed  but  only  a  3-point  or  2-point  interpolation  is  then  used. 

The  UMPE  model  also  allows  the  user  to  request  the  transmission  loss  to  the  surface.  We 
choose  to  define  this  TL  value  differently  than  Eq.  (9.1)  since  obviously  this  is  undefined  at 
(z  -  0)  -  0.  Its  derivation  is  analogous  to  the  calculation  of  rough  surface  scatter  in  section 
4.2. 1  where  it  was  shown  that  the  first  non-vanishing  moment  of  the  field  is  given  by 


<Vm>- 

*o 


Eft*1*-*) 


(9.2) 


where  we  have  included  the  scaling  factor  1  /k\  to  produce  a  quantity  with  dimensions  of  tj>2 .  We 
therefore  define  the  surface  transmission  loss  in  terms  of  the  dimensionless  vertical  derivative  of 
the  field,  i.e. 


l|.5(r)  -  (9.3) 

This  is  effectively  the  measure  of  the  field  observed  by  a  dipole  receiver.  Considering  the  wave- 
number  representation  as  in  section  6,  we  see  that  (9.3)  may  be  computed  analytically  by 


where  the  summation  is  over  all  of  wavenumber  space. 


(9.4) 


10.  Future  upgrades 

In  the  preceding  sections,  the  physics  of  the  UMPE  model  has  been  derived  and  numerical 
algorithms  have  been  developed  to  produce  solutions  of  the  underwater  acoustic  pressure  field. 
The  ultimate  goal  of  any  acoustics  model  is  the  ability  to  take  only  environmental  and  source 
parameter  inputs  and  return  an  accurate  solution.  However,  as  is  typical  of  most  numerical  models, 
several  ambiguous  variables  have  been  introduced  into  the  UMPE  model.  Default  values  have 


57 


been  defined  for  all  but  the  reference  sound  speed  which  would  hopefully  lead  to  an  efficient  and 
accurate  result.  A  method  of  computing  a  default  reference  sound  speed  is  described  below  which 
would  eliminate  the  need  for  a  “best  guess”  by  the  user.  The  corresponding  accuracy  of  the  solu¬ 
tion  provided  by  this  method  is  still  unclear. 

Other  upgrades  to  consider  might  follow  from  the  introduction  of  new  physics  (e  g.,  inter¬ 
nal  wave  scattering,  effects  of  currents,  backscattering  or  reverberation).  The  ability  to  compute 
die  full  4-D  field  (3  spatial  dimensions  and  time)  is  currently  being  investigated.  To  achieve  this 
requires  the  inclusion  of  azimuthal  coupling  in  the  UMPE  model.  A  basic  description  of  the  for¬ 
mulation  of  such  a  model  is  also  given  below. 

Finally,  one  may  consider  various  improvements  in  the  numerical  implementation  that 
increase  the  speed  and  efficiency  of  the  calculations.  The  most  obvious  extension  of  the  model  is 
a  vectorized  version  of  the  source  code  allowing  compilation  on  an  array  processor.  This  is  also 
currently  being  developed.  A  large  increase  in  speed  may  eventually  allow  implementation  of  the 
CQ-insensitive  model  mentioned  in  the  first  section  of  this  report.  However,  the  desire  to  compute 
long-range,  multiple  realization,  multiple  frequency,  multiple  coupled-azimuth  solutions  which 
would  tend  to  absorb  any  increase  in  speed  may  outweigh  the  desire  for  what  may  be  only  slight 
improvements  in  accuracy. 

10.1  Automatic  c0  selection 

Recall  that  the  choice  of  Cq,  hence  k$,  scales  the  wavenumbers  or  propagation  angles 
according  to  k  -  £0sin6 .  It  also  serves  to  scale  the  acoustic  index  of  refraction 
c0 

w(r,  z)  -  — -r-  and,  consequently,  the  environmental  potential  function  U j  (r,  z) .  Tappert 

c  (r,  z) 

(1991c)  has  suggested  one  way  to  define  a  default  reference  sound  speed  by  requiring 


where  zb  is  the  maximum  depth  of  the  water  column.  From  this  expression  it  is  obvious  that  the 
reference  sound  speed  may  now  be  range-dependent.  If,  for  example,  we  choose  to  employ  die 
standard  PE  definition 


c0  \ 


w-2V-7tJ 


then  Eq.  (10. 1)  leads  to  the  defining  equation  for  c0  given  by 


(10.2) 


1  _  1  f  dz 

4  2t>l  c2  (r,  z) 

If  instead  we  choose  the  wide-angle  PE  definition 


(103) 


the  defining  equation  for  c0  becomes 

** 

J_  _  1  f  dz 

c0  "  2b{c(r'2) 

Finally,  employing  the  LOGPE  definition  for  the  environmental  potential. 


UAr, 


■  A&r  j 


leads  to 


(10.4) 


(10.5) 


(10.6) 


/»[c0]  “  —  Jln[c(r,z)]dz 


(10.7) 


If  the  sound  speed  profile  or  the  bottom  bathymetry  is  range-dependent,  these  definitions 
produce  a  range-dependent  reference  sound  speed.  Consequently,  the  reference  wavenumber  is 

now  also  range-dependent  and  so  is  the  phase  space  potential  function  Top  (j-)  Thus,  this 

approach  requires  the  recalculation  of  the  phase  space  propagator  each  time  c0  is  updated,  increas¬ 
ing  the  total  run-time.  Again,  a  trade-off  between  speed  and  accuracy  may  be  required. 


10.2  Azimuthal  coupling 


To  develop  a  model  that  incorporates  azimuthal  coupling,  we  begin  by  returning  to  the  orig¬ 
inal  wave  equation  in  cylindrical  coordinates  defined  by  Eq.  (1.2)  without  a  source  term. 
Substituting  (1.31)  for  the  pressure  field  and  keeping  all  terms  leads  to 


d2u  d2u  1  d  u  T  i2  2  /  \  11  A 

-y +  -T+-2  —  +  <r’-)  +— 7  w  "  0 

dr2  dz2  r  dw~  [_  4r  J 


(10.8) 


Previously,  both  terms  with  the  factors  -z  were  dropped  in  the  far-fieid  approximation  and  azi¬ 


muthal  coupling  was  considered  insignificant.  If  azimuthal  coupling  is  important,  the 
corresponding  far-field  expression  is  simply 

d2u  d2u  1  d2u  ,2  2,  ,  A 

+  +  -_  +  k Qrt  (r,  z)  u  -  0 

dr  dz2  r2d( p“ 


As  before,  we  introduce  the  operators 


P  m  A 
op  dr 


(109) 


(10.10) 


Q'  »  \n2  +  l( JL  +  J — #1.) 
Qop  L  k2W  r2dq>2) 


so  (10.9)  can  be  written  in  the  more  general  form 


(10.11) 


(10.12) 


Factoring  (10.12)  into  incoming  and  outgoing  waves  produces  the  defining  first  order  differential 


equation  for  the  outgoing  energy 


-ikn  —  -  Q'u 
0  dr  ^  °P 


(10.13) 


We  again  define 


e  -  n  -  1 


(10.14) 
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(10.15) 


and  introduce 


such  that 


1  d2 


(10.16) 


Q'op  “  (|*  +  e  +  v  +  1) 1/2  .  (10.17) 

If  we  assume  that  the  effect  of  azimuthal  coupling  is  relatively  small  (since 

■fc  (r,  z,  <p)  « -§-c  (r,  s,  cp) ,  in  general)  then  we  may  expect  v  «  p  +  e.  Therefore,  (10.17)  may 
dip  02 

be  approximated  by 

Q'„p  -  +  +  -  e„,+  iv  .  (10.18) 

We  now  define  the  PE  field  function  ip  (r,  z,  «p)  according  to  Eq.  (1 .45)  and  obtain 

Tr  -  '*•<&,♦  jy- *>♦  -  -ik°{T°r*  (,019) 

where  the  last  step  follows  by  replacing  Qop  by  the  kinetic  and  potential  energy  operators  as  in  sec¬ 
tion  1 .  These  operators  can  now  take  any  of  the  approximate  forms  previously  defined. 

The  effect  of  including  azimuthal  coupling  is  to  add  a  new  operator  to  the  propagation  equa¬ 
tion  which  we  shall  define  as 


(10.20) 


Analogous  to  the  operator  Top  in  r-space,  this  operator  is  not  diagonal  in  tp-space  and  is  therefore 
not  a  simple  multiplication  operator.  Just  as  Top  coupled  vertical  angles  of  propagation,  pro¬ 
duces  a  coupling  of  azimuthal  angles.  Following  foe  same  reasoning  as  before,  we  define  a  variable 
5  as  the  “wavenumber”  analog  to  the  azimuthal  bearing  <p  according  to 

ip(r,z,(p)  -  fty(ryz,  k0s)eko3i(,d(k0s)  -  FFT  [^(r.z,*^)]  .  (10.21) 
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mmjM  if 

The  propagator  that  couples  azimuths  is  then  e  0  *  ’  where 

2 

l\p(r,s)  -  (10  22) 

2  r 

As  before,  we  assume  the  commutator  cf  Uop  and  is  negligible.  Note,  however,  that  the  com¬ 
mutator  of  Top  and  Vop  is  identically  zero  so  the  order  of  their  operations  is  arbitrary. 

We  now  must  examine  more  closely  the  explicit  calculation  of  the  variables.  To  do  so,  we 
continue  with  the  analogy  to  the  vertical  wavenumber  variable  A.  The  incremental  values  of  verti¬ 
cal  wavenumber  were  defined  by 

A-wAA,  (10.23) 

where  AT  is  the  transform  size  of  the  FFT  between  z-space  and  A-space,  and  with 

AA  -  -  -  £  (10.24) 

zT  H 

where  ~j  -  2H  is  foe  total  computational  depth.  (The  depth  of  the  water  column  plus  bottom  lay¬ 
ers  is  H  and  the  total  computational  depth  includes  the  real  and  image  ocean  or  2H.)  Since 
2H 

Az  -  -rr,  it  follows  that 

JY 

AAA:  -  ^  .  (10.25) 

N 

Therefore,  if  we  define  M  as  the  transform  size  between  qp-spacc  and  s-space  then 

s  -  mAs,  m  -  -(^-1),^  (10.26) 

with 

*  As  -  —  (10.27) 

0  qpT 

and  qpT  is  the  total  azimuthal  interval  in  radians  (e.g.,  2n,  * ,  etc.).  The  incremental  values  of  azi¬ 
muth  are  obviously  given  by 
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so  the  azimuthal  analog  to  (10.25)  is 


<PT 

N 


(10.28) 


Aq)  - 


*oAsAq>  -  —  (10.29) 

To  compute  the  3-D  propagation  of  a  point  source,  we  simply  create  identical  starting  fields 
for  each  bearing  at  qp  -  mAqp .  f  or  q>T  <  2x,  the  entire  z-space  starting  fields  should  be  tapered 

M  _ 

near  the  ends,  or  outermost  bearings,  near  qp  -  *— Aqp.  A  2-D  FPT  between  (r,  z,  qp)  and 

(r,  k,  s)  is  then  used  in  conjunction  with  a  split-step  algorithm  analogous  to  Eq.  (1.29)  with  the 
operator  substitutions 

(/op(r,z)-C/op(r,z,q>)  (10.30) 

and 

Top(r,k)^Top(r,k)  +Vop(r,s)  (10.31) 

The  above  discussion  outlines  the  basic  formulation  of  a  fully  3-D  PE/SSF  model.  Natu¬ 
rally,  there  are  subtle  issues  that  complicate  the  implementation  of  these  equations.  These  are 
currently  being  addressed  and  it  is  hoped  that  a  successful  source  code  will  be  available  within  the 
next  year. 

II.  Numerical  implementation  and  organization 

The  source  code  for  the  UMPE  model  has  been  separated  into  various  subroutines  which 
are  tasked  with  the  major  components  of  the  SSF  algorithm.  The  names  of  the  main  program  and 
subroutines  with  brief  descriptions  are  given  below.  Each  has  standard  Fortran  77  form  and  can  be 
compiled  on  Sun  workstations  using  the  included  “makefile”  script. 

pemp.f  -  This  is  the  main  program  of  the  UMPE  model.  It  handles  the  majority  of  the 
input/output.  It  also  forms  the  basic  SSF  algorithm  and  implements  this  over  many  range  steps  and 
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for  multiple  frequencies  and  multiple  rough  interface  realizations  if  required.  The  surface  interface 
boundary  condition  is  also  imposed  within  this  program 

initk.f  -  The  source  function  is  created  in  the  wavenumber  domain.  The  various  source 
parameters  required  for  this  calculation  are  read  by  the  main  program  and  passed  to  this  subroutine 
through  a  common  block. 

phase  I  f  -  The  k-domain  propagator  function  is  computed  and  passed  to  the  main  program . 
The  tapering  filter  applied  over  the  outer  1/8  of  the  positive  and  negative  wavenumbers,  and  the 
real  and  image  ocean  depths,  is  created  in  this  subroutine  and  used  in  the  creation  of  the  z-domam 
propagator.  Because  this  version  of  the  UMPE  model  assumes  constant  range  step  sizes  and  a  con¬ 
stant  reference  sound  speed  c0,  this  function  does  not  change  with  range  and  can  be  computed  prior 
to  entering  the  main  range  loop. 

phase2.f  -  The  z-domain  propagator  function  is  computed  and  passed  to  the  main  program. 
Because  this  function  contains  all  the  information  about  the  environment,  the  environmental  data 
is  input  exclusively  to  this  subroutine.  This  reduces  the  need  to  pass  data  into  the  subroutine  from 
the  main  program.  The  depth  interpolation  of  the  sound  speed  profile  as  well  as  the  range  interpo¬ 
lation  of  the  environment  is  performed  here.  If  shear  is  present  in  either  bottom  layer,  equivalent 
fluid  properties  are  computed.  The  water  column,  sediment  layer,  and  basement  layer  are  then 
combined  to  form  the  environmental  potential  function  and,  subsequently,  the  z-spacs  propagator. 

ssi.f  -  Called  by  phase2.f,  this  subroutine  performs  the  depth  interpolation  of  the  sound 
speed  profile  using  a  simple  1-2-1  filter.  This  smoothing  may  be  enhanced  by  increasing  the  num¬ 
ber  of  iterations  of  application  of  the  filter. 

etagen.f  -  Called  by  phase2  .f  if  stochastically  rougn  interfaces  are  requested.  For  exact 
surface  scatter,  the  first  derivative  of  the  surface  roughness  is  also  computed. 

ftime.f  -  When  a  broadband  analysis  is  requested,  the  main  program  outputs  to  a  single  file 
the  complex  field  function  ip  (r)  at  the  appropriate  range  for  each  frequency.  This  subroutine  is 
called  by  the  main  program  after  all  frequencies  have  been  computed  to  perform  the  necessary  Fou¬ 
rier  analysis  and  output  the  travel  time  calculations. 
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fft842.f  -  Performs  ail  of  the  Fourier  transforms  required  in  the  model.  It  computes  only 
a  forward  FFT  The  maximum  complex  array  size  allowed  is  4096  points  as  written  but  can  be 
upgraded  to  laiger  array  sizes  if  one  is  careful  to  be  consistent  throughout  the  source  code. 

tloss.f  -  Determines  values  of  TL  at  arbitrary  depths  by  performing  a  five-point  interpola¬ 
tion  of  the  field.  The  square  modulus  of  the  field  function  (r)  at  any  given  depth  is  passed  back 
to  the  main  program. 

ave.f  -  Performs  t  ange  averages  of  several  output  data  options  at  the  end  of  the  main  pro¬ 
gram. 

datime.f,  qtime.f,  qdate.f  -  Used  for  timing  purposes,  these  subroutines  allow  the  main  pro¬ 
gram  to  output  a  time  and  date  stamp  at  the  beginning  and  end  of  each  run. 

There  are  three  input  files  for  the  UMPE  model.  The  names  of  these  files  with  brief  descrip¬ 
tions  are  given  below.  Sample  listings  can  be  found  in  the  following  section  of  examples.  Only 
the  name  of  the  first  file  must  remain  the  same  as  it  is  assumed  by  the  main  program. 

pefiles.dat  -  This  file  contains  the  names  of  the  various  I/O  files  and  flags  indicating  which 
data  should  be  computed  and  output.  The  first  two  files  named  contain  the  various  run  parameters 
for  the  model  and  the  environmental  data,  respectively  (these  names  may  be  changed).  These  will 
be  discussed  in  more  detail  below.  The  remaining  filenames  define  the  names  of  various  output 
files  to  be  created  by  the  main  program  (these  names  may  also  be  changed).  Each  filename  or  set 
of  associated  filenames  is  preceded  by  one  or  more  integers.  The  first  integer  indicates  whether  the 
following  data  file  should  be  computed  and  output  (0=no,  l=yes).  The  second  integer,  if  it  exists, 
specifies  the  frequency  number  (during  broadband  calculations)  for  which  this  calculation  should 
be  made.  The  third  and  fourth  integers,  if  either  exists,  specifies  the  roughness  realization  and 
range  step  number,  respectively,  for  which  this  calculation  should  be  made.  A  set  of  filenames  are 
each  associated  with  the  same  function  but  are  defined  to  display  different  characteristics.  In  the 
order  that  they  appear  in  pefiles.dat,  the  output  filenames  and  a  brief  description  of  their  content 
are  as  follows. 
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tlb.dat,  tlbave.dat  -  TL  to  the  bottom  (basement)  interlace  and  its  associated  range 
average.  The  output  is  a  two  column  ascii  file  of  range,  TL  values.  The  range  units 
are  [km]  while  the  TL  units  are  [dB  re  lm] 

tlbs.dat,  tlbsave.dat  -  TL  to  the  bottom  sediment  interface  and  its  associated  range 
average.  The  output  is  a  two  column  ascii  file  of  range,  TL  values.  The  range  units 
are  [km]  while  the  TL  units  are  [dB  re  lm]. 

tls.dat,  tlsave.dat  -  TL  to  the  surface  interface  and  its  associated  range  average 
(as  defined  by  Eqs.  (9.1)  and  (9.4)).  The  output  is  a  two  column  ascii  file  of  range, 
TL  values.  The  range  units  are  [km]  while  the  TL  units  are  [dB  re  lm]. 
apvelrr.dat  -  The  real  part  of  the  range  component  of  acoustic  particle  velocity. 
The  output  is  a  three  column  ascii  file  of  field  variables  range,  depth,  and  acoustic 
particle  velocity.  The  range  units  are  [km],  the  depth  units  are  [m],  and  the 
velocity  units  are  [m/s]. 

apvelri.dat  -  The  imaginary  part  of  die  range  component  of  acoustic  particle 
velocity.  The  output  is  a  three  column  ascii  file  of  field  variables  range,  depth,  and 
acoustic  particle  velocity.  The  range  units  are  [km],  the  depth  units  are  [m],  and  the 
velocity  units  are  [m/s]. 

apvelzr.dat  -  The  real  part  of  the  depth  component  of  acoustic  particle  velocity. 
The  output  is  a  three  column  ascii  file  of  field  variables  range,  depth,  and  acoustic 
particle  velocity.  The  range  units  are  [km],  the  depth  units  are  [m],  and  the 
velocity  units  are  [m/s]. 

apvelzi.dat  -  The  imaginary  part  of  the  depth  component  of  acoustic  particle 
velocity.  The  output  is  a  three  column  ascii  file  of  field  variables  range,  depth,  and 
acoustic  particle  velocity.  The  range  units  are  [km],  the  depth  units  are  [m],  and  die 
velocity  units  are  [m/s]. 

prssrr.dat  -  The  real  part  of  the  complex  acoustic  pressure.  The  output  is  a 
three  column  ascii  file  of  field  variables  range,  depth,  and  acoustic  pressure.  The 
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range  units  are  [km],  Ae  depth  units  are  [m],  and  the  pressure  units  are  [Pa]. 
prssri.dat  -  The  imaginary  part  of  the  complex  acoustic  pressure.  The  output  is  a 
three  column  ascii  hie  of  field  variables  range,  depth,  and  acoustic  pressure  The 
range  units  are  [km],  die  depth  units  are  [m],  and  the  pressure  units  are  [Pa]. 
phaselr.dat,  phaseli.dat  -  The  real  and  imaginary  parts  of  the  k-domain  propagator 
function.  The  output  is  a  two  column  ascii  file  of  wavenumber,  function  values. 
The  wavenumber  units  are  [1/m].  (The  propagator  is  unitless.) 
phase2r.dat,  phase2i.dat,  ssp.dat,  potent.dat  -  The  real  and  imaginary  parts  of  the 
z-domain  propagator,  the  total  interpolated  sound  speed  profile,  and  the  total 
potential  function.  Each  output  is  a  two  column  ascii  file  of  depth,  function  values. 
The  depth  units  are  [m],  the  sound  speed  units  are  [m/s].  (The  propagator  and 
potential  are  unitless.) 

psikr.dat,  psiki.dat,  psizr.dat,  psizi.dat,  psiztdat  -  The  real,  imaginary  parts  of  the 
k-domain  field  function  i]>  (k)  and  the  real,  imaginary,  and  total  amplitude  of  the 
z-domain  field  function  ip  (z) .  Each  output  is  a  two  column  ascii  file  of  either 
wavenumber  or  depth,  function  values. 

tlgmt.dat,  botgmt.dat,  sedgmt  dat,  surfgmt.dat  -  The  complete  range-depth  field 
of  TL  data,  the  bottom  (basement)  interface  depth,  the  sediment  interface  depth, 
and  the  surface  interface  depth.  The  first  hie  is  a  three  column  ascii  file  of  field 
variables  range,  depth  and  TL.  The  remaining  files  are  two  column  ascii  files 
of  range,  interface  depth  values. 

finfield.dat  -  A  two  column  file  specifying  the  real  and  imaginary  parts  of  the 
field  function  ip  (r)  at  the  end  of  the  calculation. 


The  first  two  files  named  in  pefiles.dat  contain  all  the  input  data  for  the  UMPE  model.  As 
examples,  we  shall  refer  to  them  here  as  perun.dat  and  peenv.dat.  Brief  descriptions  of  their  con¬ 
tents  are  as  follows. 
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perundat  -  Note  that  each  line  has  a  brief  description  included  to  remind  the  user  what 
each  parameter  defines.  These  descriptions  should  not  be  changed.  The  code  assumes  the  data  is 
given  in  this  order  and  follows  the  33rd  column  of  each  line  (a  colon  exists  at  the  33rd  column). 
Otherwise,  the  data  is  given  in  free  format  with  the  exception  that  integer  inputs  are  assumed  to  be 
given  as  integers.  The  parameters  are: 

Title  -  An  ascii  title  to  simply  remind  the  user  of  the  problem  for  which  these 
parameters  were  used. 

Type  of  PE  approx.  (1-5)  -  The  integer  choice  corresponds  to  the  types  of 
approximations  described  in  section  1  of  this  report. 

iype  of  density  smoothing  (1,2)  -  The  integer  choice  defines  the  type  of  density 
mixing  function  to  be  employed,  either  a  hyperbolic  tangent  (1)  or  a  cubic  spline 
(2),  as  described  in  section  2.1.2. 

iype  of  starting  field  -  The  integer  choice  defines  the  type  of  starting  field  to  be 
employed,  either  a  Gaussian  (1),  a  wide-angle  (2),  or  a  vertical  array  (3),  as 
described  in  section  7. 

iype  of  input  data  units  -  The  integer  choice  defines  the  type  of  input  data  units, 
either  standard  MKS  (1)  (with  ranges  given  in  [km]  rather  than  [m])  or  English 
units  (2)  (with  ranges  given  in  [nm]  rather  than  [ft]).  Either  set  of  units  must  be 
used  consistently  throughout  both  input  data  files.  All  output  files  are  given  in 
MKS  units. 

Apply  surface  filter?  -  The  integer  choice  (0=no,  l=yes)  determines  whether  a 
surface  filter  will  be  applied  to  remove  surface  reflections  as  described  in  section  8. 
Approx,  or  exact  surf?  -  The  integer  choice  (0=approx,  Inexact)  determines 
whether  the  approximate  or  exact  method  of  computing  surface  scatter  should  be 
used  as  described  in  section  4.2. 

#  of  roughness  realizations  -  Integer  number  of  stochastic  interface  roughness 
realizations  for  which  calculations  should  be  made.  Outputs  of  TL  at  single  depths 
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or  at  interfaces  are  then  averaged  incoherently  to  yield  a  measure  of  the  mean 
square  pressure. 

Range  average  TL  curves?  -  The  integer  choice  (0=no,  l=yes)  determines  whether 
the  output  of  TL  at  single  depths  or  at  interfaces  should  be  range  averaged  and 
output  to  a  separate  file. 

Range  to  average  over  -  Real  number  defining  the  range  interval  to  perform  a 
running  range  average  if  requested  by  preceding  choice. 

SL  -  Real  source  level  in  dB. 

Source  depth  -  Real  depth  of  source. 

Ref.  sound  speed  -  Real  reference  sound  speed  to  be  used  in  calculations. 

Center  frequency  -  Real  center  frequency  of  source. 

Number  of  frequencies  -  Integer  number  of  frequencies  for  broadband 
calculations. 

Frequency  bandwidth  -  Real  frequency  bandwidth  of  source  for  broadband 
calculations. 

FFT  transform  size  -  Integer  number  defining  the  size  of  the  arrays  in  the 
PE/SSF  algorithm  (default  =  0). 

Range  step  size  -  Real  size  of  range  mesh  for  calculations  (default  =  0.0). 
Maximum  range  -  Real  max  range  of  calculation. 

Max  computational  depth  -  Real  max  depth  of  calculation. 

Central  source  angle  -  Real  angle  describing  steering  of  source  function. 

Source  array  length  -  Real  length  of  vertical  array  (used  only  when  vertical 
array  source  is  chosen). 

Rms  source  width  -  Real  angle  defining  width  of  source  function. 

Source  taper  width  -  Real  angle  defining  width  of  source  function  taper. 

Density  mix  length  -  Real  number  of  depth  meshes  defining  size  of  density 
mixing  function  (default  =  0.0). 
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Sspeed  mix  length  -  Real  number  of  depth  meshes  defining  size  of  sound  speed 
mixing  function  (default  =  0.0). 

Surf  loss  fctn  width  -  Real  number  of  depth  meshes  defining  size  of  surface 
loss  function  used  for  the  approximate  rough  surface  scatter  (default  -  0.0). 

#  TL  depths,  TL  depths  -  Integer  number  of  depths  (<=  10),  and  real  depths  that 
follow,  for  which  calculations  of  TL  should  be  computed  and  output. 

#  BB  ranges,  BB  ranges  -  Integer  number  of  ranges  (<=  4),  and  real  ranges  that 
follow,  for  which  broadband  calculations  of  time-domain  TL  field  data  should  be 
computed  and  output. 

BB  extraction  depth  -  Real  single  depth  of  TL  -vs-  time  data  to  be  extracted 
from  each  of  the  above  BBTL  ranges. 

peenv.dat  -  The  following  environmental  parameters  are  contained  in  this  file. 

ISEED  -  An  integer  seed  for  random  number  generators. 

NSS  -  Integer  number  of  different  sound  speed  profiles  to  follow  (must  be  at  least 
one).  Each  profile  has  the  following  format. 

ZV,  NSSDA  -  Real  range  of  current  profile  (first  profile  must  be  at  range 
0.0)  and  integer  number  of  sound  speed  values  in  depth  given  for  current 
profile  which  follow. 

SSD,  SS  -  Real  depth  of  this  sound  speed  value  and  real  sound 
speed  value. 

NB  -  Integer  number  of  bottom  profiles  to  follow  (must  be  at  least  one).  Each 
profile  has  the  following  format. 

ZB,  BD,  BV,  BG,  RDEN,  BLKMI,  SIGBOT,  ALCBO,  BSWS,  BSWLKMI  - 
Range  of  this  profile,  depth,  sound  speed  at  top  of  this  profile,  sound  speed 
gradient  for  this  profile,  density  ratio  w.r.t.  water,  attenuation,  rms  roughness 
of  interface,  correlation  length  of  interface  roughness,  shear  wave  sound. 
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speed,  and  shear  wave  attenuation.  (All  parameters  are  real.) 

NBS  -  Integer  number  of  sediment  profiles  to  follow  (may  be  zero).  Each 
profile  has  the  following  format. 

ZBS,  BSD,  BSV,  BSG,  RDENS,  BSLKMI,  SIGBOTS,  ALCBOS,  BSSWS, 
BSSWLKMI  - 

Range  of  this  profile,  sediment  thickness,  sound  speed  at  top  of  this  profile, 
sound  speed  gradient  for  this  profile,  density  ratio  w.r.t.  water,  attenuation,  rms 
roughness  of  interface,  correlation  length  of  interface  roughness,  shear  wave 
sound  speed,  and  shear  wave  attenuation.  (All  parameters  are  real.)  Note 
that  sediment  profiles  are  defined  in  terms  of  sediment  thicknesses  and  not 
true  depth.  The  sediments  are  layered  on  top  of  the  bottom  profiles. 

NS  -  Integer  number  of  surface  roughness  profiles  to  follow  (may  be  zero).  Each 
profile  has  the  following  format. 

ZS,  SIGSUR  -  Range  of  this  profile,  wind  speed  (Internally,  wind  speed 
is  converted  into  an  rms  roughness  and  correlation  length  scale  as  described 
in  section  4.2) 

NSB  -  Integer  number  of  equivalent  bubble  surface  roughness  profiles  to  follow 
(may  be  zero).  Each  profile  has  the  following  format. 

ZSB,  SIGSURB  -  Range  and  equivalent  rms  roughness  of  this  profile. 

12.  Examples 

In  this  section,  we  present  several  examples  of  input  files  and  output  solutions.  The  first 
three  were  taken  from  problems  defined  in  the  PE  Workshop  II  proceedings  (ref?)  and  are  designed 
to  exhibit  the  accuracy  of  the  UMPE  model  when  default  values  are  selected.  Each  of  these  has 
reference  solutions  with  which  we  can  compare  the  results  of  the  model.  The  last  three  problems 
were  chosen  to  exhibit  some  of  the  versatile  features  of  the  model. 
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Example  1 :  The  first  example  taken  from  the  PE  Workshop  11  proceedings  consists  of  a 
simple  fluid  half-space.  The  solution  displays  a  simple  pattern  due  to  the  interference  of  the  unage 
source  and  the  real  source.  This  problem  is  sometimes  referred  to  as  the  Lloyd’s  mirror  problem 
It  is  listed  as  Test  Case  1  in  the  workshop  proceedings.  The  output  requested  from  the  UMPE 
model  is  the  transmission  loss  at  a  single  depth  (defined  in  perun.dat)  and  the  complete  71  field 
data  (requested  in  pefiles.dat).  Note  that  since  die  UMPE  model  assumes  at  least  one  bottom  inter¬ 
face,  we  simply  define  the  bottom  with  fluid  properties  equivalent  to  the  overlying  fluid  In  Fig.  1 , 
the  71  field  data  is  displayed  and  the  expected  interference  pattern  can  easily  be  seen.  The  com¬ 
parison  of  the  71  reference  solution  with  the  model  prediction  is  shown  in  Fig.  2.  (All  figures  have 
been  created  with  the  GMT  plotting  package  (ref).) 


Input  files  for  Test  Case  1  of  PE  Workshop  II: 


TITLE  :PEH  lest  case  1 

TYPE  OF  PE  APPROX.  (1-5)  :4 

TYPE  DENSITY  SMOOTHING  (U)  :2 

TYPE  OF  STARTING  FIELD  (1-3)  :2 

TYPE  OF  INPUT  DATA  UNITS  (1,2)  :  1 

APPLY  SURFACE  FILTER?  (0,1)  :0 

APPROX  OR  EXACT  SURF?  (0,1 )  :0 

#  OF  ROUGHNESS  REALIZATIONS  :  1 

RANGE  AVERAGE  TL  CURVES?  (0,1 )  :0 

RANGE  TO  AVERAGE  OVER  (km)  :  1 .0 

SL(dB  re  uPa re  lmjyd)  :20u.000 

SOURCE  DEPTH  (m,ft)  :350.000 

REF.  SOUND  SPEED  (m/s,ft/s)  :  1 500.00 

CENTER  FREQUENCY  (Hz)  :40.000 

NUMBER  OF  FREQUENCIES  :  1 

FREQUENCY  BANDWIDTH  (Hz)  : 32. 0000 

FFT  TRANSFORM  SIZE  :0 

RANGE  STEP  SIZE  (kmjim)  :0.0 

MAXIMUM  RANGE  (knunn)  :10.000 

MAX  COMPUTATIONAL  DEPTH  (m,ft)  :6000.00 

CENTRAL  SOURCE  ANGLE  (DEG)  :0. 

SOURCE  ARRAY  LENGTH  (m,ft)  :20.0000 

RMS  SOURCE  WIDTH  (DEG)  :80.0000 

SOURCE  TAPER  WIDTH  (DEG)  :  10.00000 

DENSITY  MIX  LENGTH  (#  DELD)  :0.0 

SSPEED  MIX  LENGTH  (#  DELD)  :0.0 

SURF  LOSS  FCTN  WIDTH  (#  DELD)  :0.0 

#TL  DEPTHS,  TL  DEPTHS  (m,ft)  :1  3990.000 

#BB  RANGES,  BB  RANGES  (knMim)  :1  100.000 
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BB  EXTRACTION  DEPTH  (m,ft) 


:  1000.000 


peenv.dat 


3000 

1 

0.02 

0.  1500.00 
10000.0  1500. 

1 

0.0  3990.0  1500.00  0.0 
1 

0.0  0.  150000  0.0  1.0 

1 

0.0  0.0 
1 

0.0  0.0 


1.0  .000  00.0  200.0  0.0  .000 
.000  00.0  200.0  0.0  .000 


Example  2:  This  problem  is  listed  as  Test  Case  2  in  the  workshop  proceedings  It  is  char¬ 
acterized  by  a  water  column  of  varying  depth  overlying  a  denser  fluid  bottom.  The  interface  forms 
an  upslope-downslope  configuration.  Only  transmission  loss  at  two  depths  is  requested  as  output. 
The  results  are  compared  with  the  reference  solutions  in  Fig.  3. 


Input  files  for  Test  Case  2  of  PE  Workshop  II: 


pcruiLdfli 


TITLE  :PED  test  ca*  2 

TYPE  OF  PE  APPROX.  (1-5)  :4 

TYPE  DENSITY  SMOOTHING  (1 2)  :2 

TYPE  OF  STARTING  FIELD  (1-3)  :2 

TYPE  OF  INPUT  DATA  UNITS  (U)  :1 

APPLY  SURFACE  FILTER?  (0,1)  :0 

APPROX  OR  EXACT  SURF?  (0,1)  :0 

*  OF  ROUGHNESS  REALIZATIONS  :  1 

RANGE  AVERAGE  TL  CURVES?  (0,1 )  :0 

RANGE  TO  AVERAGE  OVER  (km)  :  1 .0 

SL(dB  re  uPa  re  !m,lyd)  : 200.000 

SOURCE  DEPTH  (m,ft)  :  100.000 

REF.  SOUND  SPEED  (m/8,ft/s)  :  1 500.00 

CENTER  FREQUENCY  (Hz)  :25.000 

NUMBER  OF  FREQUENCIES  :  1 

FREQUENCY  BANDWIDTH  (Hz)  :  32.0000 

FFT  TRANSFORM  SIZE  :0 

RANGE  STEP  SIZE  (kmjsn)  :0.0 

MAXIMUM  RANGE  (km,nm)  :7.000 

MAX  COMPUTATIONAL  DEPTH  (m,ft)  :800.00 

CENTRAL  SOURCE  ANGLE  (DEG)  :0. 
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SOURCE  ARRAY  LENGTH  (m,tt) 

RMS  SOURCE  WIDTH  (DEG) 

SOURCE  TAPER  WIDTH  (DEG) 

DENSITY  MIX  LENGTH  (#  DELD) 

SSPEED  MIX  LENGTH  (#  DELD) 

SURF  LOSS  FCTN  WIDTH  (#  DELD) 

#  TL  DEPTHS,  TL  DEPTHS  (m,ft) 

#  BB  RANGES,  BB  RANGES  (kmjun) 

BB  EXTRACTION  DEPTH  (mjt) 

pccnv.dat 

3000 
1 

0.02 

0.  1300.00 
10000.0  1300. 

3 

0.0  200.0  1700.00  0.0  1.3  294  00.0  200.0  0.0  .000 
3.3  23.0  1700.00  0.0  1.3  .294  00.0  200.0  0.0  .000 

7.0  200.0  1700.00  0.0  1.3  294  00.0  200.0  0.0  .000 
1 

0.0  0.  1300.00  0.0  1.0  .000  00.0  200.0  0.0  .000 

1 

0.0  0.0 
1 

0.0  0.0 


:20.0000 
: 60.0000 
:20.00000 
.-0.0 
:0.0 
:0.0 

2  20.000  130.000 
:1  100.000 
:  1000.000 


Example  3:  The  problem  is  listed  as  Test  Case  7  in  the  workshop  proceedings  and  is  some¬ 
times  referred  to  as  the  “Porter’s  duct”  problem.  The  environment  consists  of  a  typical  deep  ocean 
sound  speed  channel  with  a  shallow  surface  duct.  The  bottom  is  a  simple  fast  fluid.  The  output 
requested  from  the  UMPE  model  is  the  transmission  loss  at  a  single  depth  and  the  complete  TL  field 
data.  Additionally,  we  requested  as  output  the  environmental  propagator  functions.  This  allows 
us  to  plot  the  extrapolated  sound  speed  profile  used  by  the  model.  In  Fig.  4,  the  TL  field  data  is 
displayed  and  the  sound  speed  profile  is  shown  in  Fig.  5.  The  comparison  of  the  TL  reference  solu¬ 
tion  with  the  model  prediction  is  shown  in  Fig.  6  for  two  different  values  of  reference  sound  speed. 
For  this  particular  problem,  the  solution  is  found  to  be  highly  sensitive  to  changes  in  c0  when  the 
wide  angle  approximation  is  chosen. 


Input  files  for  Test  Case  7  of  PE  Workshop  II: 
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CfiOULdftl 


TITLE  PEG  test  cmc  7 

TYPE  OF  PE  APPROX.  (1-5)  :4 

TYPE  DENSITY  SMOOTHING  (U)  :2 

TYPE  OF  STARTING  FIELD  (1-3)  :2 

TYPE  OF  INPUT  DATA  UNITS  (1,2)  :  I 

APPLY  SURFACE  FILTER?  (0,1)  :0 

APPROX  OR  EXACT  SURF?  (0,1 )  :0 

#  OF  ROUGHNESS  REALIZATIONS  :  1 

RANGE  AVERAGE  TL  CURVES?  (0,1)  .0 

RANGE  TO  AVERAGE  OVER  (km)  :1 .0 

SL(dB  re  uPn  re  lm,lyd)  : 200.000 

SOURCE  DEPTH  (nMl)  :25.000 

REF.  SOUND  SPEED  (m/s,fl/8)  :  1483.00 

CENTER  FREQUENCY  (Hz)  : 80.000 

NUMBER  OF  FREQUENCIES  :1 

FREQUENCY  BANDWIDTH  (Hz)  : 32.0000 

FFT  TRANSFORM  SIZE  .0 

RANGE  STEP  SIZE  (knyun)  :0.0 

MAXIMUM  RANGE  (kn^mn)  :  150.000 

MAX  COMPUTATIONAL  DEPTH  (m,ft)  : 6000.00 

CENTRAL  SOURCE  ANGLE  PEG)  :0. 

SOURCE  ARRAY  LENGTH  (m,ft)  :20.0000 

RMS  SOURCE  WIDTH  PEG)  : 60.0000 

SOURCE  TAPER  WIDTH  (MG)  20.00000 

DENSITY  MIX  LENGTH  (#  DELD)  :0.0 

SSPEED  MIX  LENGTH  (#  DELD)  :0.0 

SURF  LOSS  FCTN  WIDTH  (#  DELD)  0.0 

#  TL  DEPTHS,  TL  DEPTHS  (m,ft)  :  1  100.000 

#  BB  RANGES,  BB  RANGES  (kmjim)  :  1  100.000 

BB  EXTRACTION  DEPTH  (m,ft)  :  1000.000 


peenv.dat 

3000 

1 

0.0  20 
0.  1497.00 


250. 

1502.00 

300. 

1485.00 

375. 

1478.00 

425. 

1477.00 

500. 

1476.00 

600. 

1476.50 

700. 

1477.00 

810. 

1478.00 

900. 

1479.00 

1000. 

1480.00 

1100. 

1481.00 

1180. 

1482.00 

1340. 

1484.00 

1600. 

1487.00 

1800. 

1490.00 

2500. 

1498.70 

3000. 

1506.80 

4000. 

1523.90 

75 


10000.0  1523-90 
1 

0.0  4000.0  1523.80  0.0  1.0  .066  00.0  200.0  0.0  .000 
1 

0.0  0.  1500.00  0.0  1.0  .000  00.0  200.0  0.0  .000 

1 

0.0  0.0 

1 

0.0  0.0 


Example  4:  This  example  is  chosen  to  display  the  effects  of  rough  bottom  interfaces  on 
acoustic  propagation.  We  have  defined  the  environment  as  a  typical  deep  ocean  water  column  with 
a  thin  sediment  layer  overlying  a  hard  bottom.  The  average  sediment  depth  is  10  m.  The  rms 
roughness  of  the  sediment  interface  is  5  m  while  the  rms  roughness  of  the  bottom  interface  is  20 
m.  Both  roughness  spectra  have  been  filtered  from  10  km  to  20  km  to  remove  length  scales  longer 
than  20  km  (this  number  is  set  within  the  subroutine  etagen.f).  The  TL  field  data  has  been  output 
and  is  plotted  in  Fig.  7.  The  transmission  loss  to  each  interface  has  been  averaged  over  1  km  inter¬ 
vals  and  these  curves  are  displayed  in  Fig.  8. 


Input  files  for  Example  4: 


TITLE 

TYPE  OF  PE  APPROX.  (1-5) 

TYPE  DENSITY  SMOOTHING  (12) 
TYPE  OF  STARTING  FIELD  (1-3) 

TYPE  OF  INPUT  DATA  UNITS  (12) 
APPLY  SURFACE  FILTER?  (0,1) 
APPROX  OR  EXACT  SURF?  (0,1) 

#  OF  ROUGHNESS  REALIZATIONS 
RANGE  AVERAGE  TL  CURVES?  (0,1) 
RANGE  TO  AVERAGE  OVER  (km) 
SL(dB  re  uPa  re  lm,lyd) 

SOURCE  DEPTH  (mjt) 

REF.  SOUND  SPEED  (m/sjt'i) 

CENTER  FREQUENCY  (Hz) 

NUMBER  OF  FREQUENCIES 
FREQUENCY  BANDWIDTH  (Hz) 

FFT  TRANSFORM  SIZE 
RANGE  STEP  SIZE  (kmjun) 
MAXIMUM  RANGE  (knunn) 

MAX  COMPUTATIONAL  DEPTH  (m,ft) 
CENTRAL  SOURCE  ANGLE  (DEG) 
SOURCE  ARRAY  LENGTH  (m,ft) 


Rough  bottom  interfaces 
:4 
2 
2 
1 

.0 

:0 

:1 

1 

1.0 

200.000 

250.000 

1500.00 

100.000 

1 

32.0000 

1024 

:0.01 

100.000 

6000.00 

0. 

20.0000 
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RMS  SOURCE  WIDTH  (DEG)  60.0000 

SOURCE  TAPER  WIDTH  (DEG)  20.00000 

DENSITY  MIX  LENGTH  (#  DELD)  0.0 

SSPEED  MIX  LENGTH  (#DELD)  :0.0 

SURF  LOSS  FCTN  WIDTH  (#  DELD)  0.0 

#  TL  DEPTHS,  TL  DEPTHS  (m,ft)  1  100.000 

#  BB  RANGES,  BB  RANGES  (kmjun)  :  1  100.000 

BB  EXTRACTION  DEPTH  (nuft)  :  1000.000 


ccsnYial 

3000 

1 

0.0  30 

O.OOOOOOE+OO  1.520840E+03 
1  .OOOOOOE+Ol  1.520660E+03 
2. 000000E+01  1.520280E+03 
3.000000E+01  1.519330E+03 
5.000000E+01  1.512440E+03 
7.500000E+01  1  505780E+03 
1.000000E+02  1.503170E+03 
1250000E+Q2  1.501060E+03 
1.500000E+Q2  1.498310E+03 
2.000000E+Q2  1.491 WOE +03 
2.500000E+02  1489080E+03 
3.000000E+02  1.486760E+03 
4.G00000E+02  1.482140E+03 
5.000000E+Q2  1.479350E+03 
6000000E+G2  1.478510E+03 
7.000000E+02  1.478650E+03 
8.000000E+Q2  1.479040E+03 
9.000000E+02  1.479640E+03 
1.000000E+03  1  480270E+03 
1.100000E+03  1.481050E+03 
1.200000E+03  1.481890E+03 
1.300000E+03  1.482820E+03 
1  400000E+03  1.483740E+03 
1.500000E+03  1.484660E+03 
1.750000E+03  1.487460E+03 
2.000000E+03  1.490590E+03 
2.500000E+03  1.497990E+O3 
3.000000E+03  1.506030E+03 
3.500000E+03  1.514420E+03 
1.000000E+04  1.514420E+03 
1 

0.0  3500.0  2000.00  0.0  3.0  .100  50.0  1000.0  0.0  .000 
1 

0.0  500.0  1700.00  0.0  1.5  .010  10.0  200.0  0.0  .000 

1 

0.0  0.0 
1 

0.0  0.0 
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Example  S:  In  this  example,  we  exhibit  the  model’s  ability  to  compute  the  exact  surface 
scatter  from  a  large  Gaussian  displacement  of  the  surface.  (This  displacement  was  defined  analyt¬ 
ically  and  included  by  specifically  editing  parts  of  the  phase2.f  subroutine  This  is  not  available  in 
the  generic  version  but  does  illustrate  how  one  may  easily  modify  the  code  to  include  various  fea¬ 
tures.)  Note  that  this  calculation  could  not  be  performed  properly  with  the  approximate  method 
since  the  small  surface  displacement  approximation  is  invalid  The  water  column  is  treated  as  a 
homogeneous  half-space  with  the  “bottom”  having  properties  identical  to  the  overlying  fluid  In 
Fig.  9  we  display  the  TL  field  data,  and  the  surface  transmission  loss,  TL^  as  defined  by  Eqs.  (9.1) 
and  (9.4),  is  shown  in  Fig.  10. 


Input  files  for  Example  S: 


Bcmn.dat 


TITLE 

TYPE  OF  PE  APPROX.  (1-5) 

TYPE  DENSITY  SMOOTHING  (U) 
TYPE  OF  STARTING  FIELD  (1-3) 

TYPE  OF  INPUT  DATA  UNITS  (1 2) 
APPLY  SURFACE  FILTER?  (0,1) 
APPROX  OR  EXACT  SURF?  (0,1) 

U  OF  ROUGHNESS  REALIZATIONS 
RANGE  AVERAGE  TL  CURVES?  (0,1) 
RANGE  TO  AVERAGE  OVER  (km) 
SL(dB  re  uPa  re  lm.lyd) 

SOURCE  DEPTH  fm,ft) 

REF.  SOUND  SPEED  (m/s,ft/s) 

CENTER  FREQUENCY  (Hz) 

NUMBER  OF  FREQUENCIES 
FREQUENCY  BANDWIDTH  (Hz) 

FFT  TRANSFORM  SIZE 
RANGE  STEP  SIZE  (kmjun) 

MAXIMUM  RANGE  (km^m) 

MAX  COMPUTATIONAL  DEPTH  (m,ft) 
CENTRAL  SOURCE  ANGLE  (DEG) 
SOURCE  ARRAY  LENGTH  (m,ft) 

RMS  SOURCE  WIDTH  (DEG) 

SOURCE  TAPER  WIDTH  (DEG) 
DENSITY  MIX  LENGTH  (#  DELD) 
SSPEED  MIX  LENGTH  (#  DELD) 

SURF  LOSS  FCTN  WIDTH  (#  DELD) 

#  TL  DEPTHS,  TL  DEPTHS  (m,ft) 

#  BB  RANGES,  BB  RANGES  (km,nm) 
BB  EXTRACTION  DEPTH  (m,ft) 


Gaussian  surface  displacement 
4 
2 
2 
1 
0 
1 
1 
0 

1.0 

200.000 

75.000 

1500.00 

250.000 

1 

32.0000 

1024 

0.01 

5.000 

1500.00 

0. 

20.0000 

40.0000 

20.00000 

0.0 

0.0 

0.0 

0  100.000 
1  lOO.O'JO 
1000.000 
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txcny.dai 


3000 

1 

0.02 

0.0  1300.00 

10000.0  1500.00 
1 

0.0  1200.0  1500.00  0.0  1.0  .000  0,0  200.0  0.0  .000 
1 

0.0  0.0  1500.00  0.0  1.0  .000  0.0  200.0  0.0  .000 

1 

0.0  0.0 
1 

0.0  0.0 


Example  6:  In  our  final  example,  the  broadband  capabilities  of  the  model  are  displayed. 
The  environment  we  have  chosen  consists  of  a  typical  deep  ocean  sound  speed  profile  overlying  a 
hard,  lossy  bottom  and  is  range-independent.  As  output  we  have  selected  the  time  domain  TL  field 
and  a  slice  at  the  sound  channel  axis  at  600  m.  The  normalization  of  the  TL  levels  was  performed 
after  the  calculation  with  the  peak  level  after  one  range  step.  Because  a  range  step  is  10  m,  foe 

values  are  rescaled  by  subtracting  TLmin  (r  -  10  m )  -  lOlog  (r  "  J9.  m) .  The  number  of  fre- 

quencies  was  chosen  to  yield  a  time  window  width  of  roughly  5  sec.  A  more  narrow  source 
function  was  defined  to  avoid  high  angle  propagation  that  would  arrive  significantly  later  than  foe 
water-bourne  eneigy.  The  time  domain  TL  field,  displayed  in  Fig.  11,  exhibits  foe  typical  caustic 
structure  for  deep  sound  channel  propagation.  The  first  figure  shows  foe  full  time  window  and  foe 
second  expands  foe  view  to  examine  the  earliest  arrivals.  Source  angles  beyond  roughly  1 5  °  have 
long  enough  path  lengths  to  introduce  relatively  long  travel  times.  Such  arrivals  can  be  seen  to 
arrive  as  much  as  4  sec  later  than  the  horizontal  axial  arrival  and  are  observed  to  create  a  wrap¬ 
around  effect  in  the  calculation.  This  provides  a  clear  example  of  the  need  for  a  sufficient  time  win¬ 
dow  width.  The  time  arrival  pattern  observed  at  the  axis  of  foe  sound  channel  is  displayed  in  Fig. 

1 2.  Again,  an  expanded  view  of  the  earliest  arrivals  is  provided.  Significant  side-lobe  eneigy  is 
present  due  to  foe  high  ratio  of  bandwidth  to  center  frequency. 
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Input  files  for  Example  6: 
perun.dat 


TITLE 

TYPE  OF  PE  APPROX.  (1-5) 

TYPE  DENSITY  SMOOTHING  (1,2) 
TYPE  OF  STARTING  FIELD  (1-3) 

TYPE  OF  INPUT  DATA  UNITS  (1 .2) 
APPLY  SURFACE  FILTER?  (0,1) 
APPROX  OR  EXACT  SURF?  (0,1) 

#  OF  ROUGHNESS  REALIZATIONS 
RANGE  AVERAGE  TL  CURVES?  (0,1) 
RANGE  TO  AVERAGE  OVER  (km) 
SL(dB  re  uPa  re  lm,lyd) 

SOURCE  DEPTH  (m,ft) 

REF.  SOUND  SPEED  (m/s,ft/s) 

CENTER  FREQUENCY  (Hz) 

NUMBER  OF  FREQUENCIES 
FREQUENCY  BANDWIDTH  (Hz) 

FFT  TRANSFORM  SIZE 
RANGE  STEP  SIZE  (Imumi) 
MAXIMUM  RANGE  (km,mn) 

MAX  COMPUTATIONAL  DEPTH  (m,ft) 
CENTRAL  SOURCE  ANGLE  (DEG) 
SOURCE  ARRAY  LENGTH  (m,ft) 

RMS  SOURCE  WIDTH  (DEG) 

SOURCE  TAPER  WIDTH  (DEG) 
DENSITY  MIX  LENGTH  (#  DELD) 
SSPEED  MIX  LENGTH  (#  DELD) 

SURF  LOSS  FCTN  WIDTH  (#  DELD) 

#  TL  DEPTHS,  TL  DEPTHS  (m,ft) 

#  BB  RANGES,  BB  RANGES  (knMun) 
BB  EXTRACTION  DEPTH  (mjt) 


Deep  ocean,  broadband,  range-independent 
:4 

2 

2 

:1 

:0 

:0 

:1 

:0 

:1.0 

:200.000 
.-600.000 
:1478.S1 
: 50.000 
:256 

: 50.0000 

:1Q24 

:0.01 

:  100.000 

:6500.00 

:0. 

:20.0000 
-.20.0000 
:  10.00000 
:0.0 
:0.0 
:0.0 

:0  100.000 
2  0.01  100.000 
:600.000 


3000 

1 

0.0  32 

O.OOOOOOE+OO 

I.000000E+O1 

2.000000E+01 

3.000000E+01 

5.000000E+01 

7.500000E+01 

1.000000E+02 

1Z50000E+02 

1.500000E+02 

2.000000E+02 

2.500000E+02 

3.000000E+02 

4.000000E+02 

5.000000E+02 

6.000000E+02 

7.000000E+02 

8.000000E+02 


1.520840E+03 
1.520660E+03 
1.520280E+03 
1.519330E+03 
1.512440E-K)3 
1.505780E+03 
1.503170E+03 
1.501060E+03 
1.498310E+03 
1.491900E+03 
1.489080E+03 
1.486760E+03 
1.482140E+03 
1 479350E+03 
1.478510E+03 
1.478650E+03 
1.479040E+03 
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9.000000E+02  1479640E+03 
1.000000E+03  1.480270E+03 
1.100000E+03  1.481050E+03 
1.200000E+03  1.481890E+03 
1.300000E+03  1  482820E+03 
l  400000E+03  1 483740E+03 
1.500000E+03  1  484660E+03 
1. 75000® +03  1.487460E+03 
2.00000® +03  1.490590E+03 
2.500000E+03  1.497990E+03 
3.000000E+03  1.506030E+03 
3.500000E+03  I.514420E+03 
4.000000E+03  1.523210E+03 
4.500000E+03  1.532230E+03 
1.000000E+04  1.532230E+03 
1 

0.0  4500.0  1700.00  0.0  2.0 
1 

0.0 
1 

0.0 
1 

0.0 


J200  00.0  200.0  0.0  .000 


0.  1500.00  0.0  1.0  .000  00.0  200.0  0.0  .000 


0.0 


0.0 
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P€U  Test  Cate  1 


Kano*  (km) 


I  Af  17  Owl  UMtLfiMNI 


Figure  1 :  Transmission  loss  field  plot  calculated  from  UMPE  model  for 
Test  Case  1  of  PE  Workshop  II. 


PEH  Tad  Case  1 


Figure  d:  Comparison  of  TL  line  plots  from  UMPE  mod  *'  (solid  curve) 
and  reference  solution  (dashed  curve)  for  Test  Case  1  of  PE  Workshop  II 


pen  Twt  cmc  2b 


Figure  3:  Comparison  of  TL  line  plots  at  two  different  depth  from  UMPE 
model  (solid  curves)  and  reference  solutions  (dashed  curves)  for  Test  Case 
2  of  PE  Workshop  H 
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TrmM«on  La**  (40  r*  \m) 

PEN  TMt  C*M  7 


Figure  4:  Transmission  loss  field  plot  calculated  from  UMPE  model  for 
Test  Case  7  of  PE  Workshop  II. 


PEII  Te«t  Com  7  Sound  Speed  Profile 

loMStMUM) 

1*»  H«  1«W  tlH  tilt 


Figure  5:  Extrapolated  sound  speed  profile  used  by  the  UMPE  model 
for  Test  Case  7  of  PE  Workshop  D. 
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ROTftOan) 

Figure  6:  Comparison  of  TL  line  plots  at  two  different  reference  sound 
speeds  from  UMPE  model  (solid  curves)  and  reference  solution  (dashed 
curves)  for  Test  Case  7  of  PE  Workshop  II. 
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Tii»  ii  nn  La  (a*  x  la) 

Rocgn  Soasrrwt  Rough  Dmwm  irtarfaoo 


Figure  7:  Transmission  loss  field  plot  displaying  effects  of  forward  scatter 
due  to  two  rough  bottom  interfaces 


Figure  8:  Plots  of  transmission  loss  computed  to  each  rough  interface. 
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Emu  Surt»e»  Samr  Guanari  StrinOwlnrani 


Figure  9:  Transmission  loss  field  plot  displaying  effects  of  exact  surface 
scattering  from  a  Gaussian  surface  displacement 


Figure  10:  Plot  of  TLS  in  presence  of  Gaussian  surface  displacement 
as  defined  by  Eqs.  (9.1)  and  (9.3)  of  the  text. 


87 


Rang*  ino«p*ooent  D— p  Ocm n  Protae 


Figure  1 1 :  Transmission  loss  field  plots  in  the  time  domain  computed 
from  a  broadband  UMPE  model  run  in  a  typical  deep  ocean  profile  The 
upper  figure  displays  the  total  time  window  computed  while  the  lower 
figure  shows  an  expanded  view  of  the  earliest  arrivals. 


4 M 

nn^jcmd  T1 mm  (— c) 


Figure  12:  Time  arrival  patterns  of  TL  at  the  channel  axis  for  a  typical 
deep  ocean  profile.  The  upper  figure  displays  die  total  time  window 
computed  while  the  lower  ngure  shows  an  expanded  view  of  the  earliest 
arrivals. 
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